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Abstract 


For  Navy  relevant  geometries,  separation  of  wall  bounded  flows  is  a  highly  complex 
phenomenon.  Because  of  the  relatively  high  Reynolds  numbers  involved,  separation  is  always 
associated  with  considerable  unsteadiness.  This  unsteadiness  is  caused  by  large  coherent 
structures  that  are  a  consequence  of  hydrodynamic  instability  mechanisms  of  the  mean  flow.  In 
addition,  due  to  the  shape  of  underwater  vehicles  (submarines,  torpedoes,  low  aspect  ratio  lifting 
or  control  surfaces)  the  separation  is  three-dimensional  (3-D).  The  combination  of  three- 
dimensionality  and  unsteadiness  results  in  a  highly  complex  time-dependent  topology  of  the 
separated  region.  In  a  combined  numerical/experimental  effort  we  investigated  laminar 
separation  bubbles  in  external  flows.  For  the  simulations  we  employed  highly-resolved  direct 
numerical  simulations  (DNS)  to  obtain  a  deeper  understanding  of  the  various  physical 
mechanism  governing  separation,  transition,  and  reattachment  of  3-D  bubbles.  Ultimately,  such 
understanding  may  pave  the  way  for  the  development  of  effective  and  efficient  strategies  for 
preventing  separation  for  practical  applications.  We  also  evaluated  hybrid  turbulence  models  for 
high  Reynolds  number  flows.  In  particular,  we  carried  out  DNS,  RANS,  and  hybrid  simulations 
of  a  turbulent  square-duct  flow.  Based  on  these  simulations  we  decided  on  two  hybrid  strategies 
for  simulating  the  asymmetric  diffuser  experiments  that  were  conducted  at  Stanford  University 
by  J.  Eaton  and  co-workers. 
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1.  INTRODUCTION 


Separation  for  Navy  relevant  geometries  (submarines,  torpedoes,  fins,  low  aspect  ratio  lifting 
or  control  surfaces)  is  often  three-dimensional  (3-D)  and  associated  with  considerable 
unsteadiness.  Strong  unsteadiness  can  be  introduced  by  coherent  flow  structures  which  can  be 
characterized  as  areas  of  organized  fluid  motion  within  a  turbulent  flow.  These  structures  are  a 
consequence  of  hydrodynamic  instabilities  and  originate  in  the  separated  shear  layer.  They 
impact  the  separation  and  reattachment  behavior  and  thereby  affect  size  and  shape  of  the 
separated  region  or  separation  bubble.  This  can  cause  considerable  unsteady  hydrodynamic  loads 
and  ultimately  a  change  in  the  lift  and  drag  characteristics  of  the  vehicle.  An  improved 
understanding  of  the  flow  physics  governing  3-D  separation  in  general  and  the  dynamics  of  such 
structures  in  particular  is  desirable  as  it  may  lead  to  techniques  that  allow  for  the  prediction  of 
unsteady  loads  which  would  in  turn  enable  the  development  of  geometries  with  more  favorable 
hydrodynamic  properties.  An  improved  understanding  may  also  point  out  possibilities  for  novel 
devices  or  strategies  that  may  help  prevent  or  control  separation  by  influencing  the  unsteady  flow 
structures. 

The  topology  of  3-D  separation  bubbles  can  be  analyzed  by  considering  the  limiting  skin 
friction  lines  [Tobak  and  Peake,  1982],  Saddle  points  are  points  where  two  skin  friction  lines 
cross  each  other.  Points  to  which  an  infinite  number  of  skin  friction  lines  converge  are  known  as 
nodes  and  foci.  Similar  structures  can  be  identified  in  wall  normal  planes,  such  as  the  symmetry 
plane.  Three-dimensional  flow  separation  cannot  be  universally  characterized  by  the  vanishing  of 
wall  shear.  Rather,  as  Lighthill  [1963]  pointed  out,  it  is  characterized  by  the  convergence  of  the 
limiting  skin  friction  lines  onto  one  separation  skin  friction  line.  The  dividing  stream-surface 
originating  from  this  line  is  called  the  separation  stream-surface.  Focal  points,  around  which  skin 
friction  lines  spiral,  are  the  roots  of  wall-normal  vortices,  often  referred  to  as  “horn  vortices 

Because  of  the  3-D  nature  of  the  flow  the  underly  ing  phy  sics  are  very  complex.  For  separated 
flows  hydrodynamic  instability  mechanisms  of  the  mean  flow  can  lead  to  coherent  structures. 
The  instability  mechanisms  may  be  3-D  and  the  resultant  structures  and  their  interaction  can  be 
very  complicated.  Progress  in  modeling  and  understanding  of  the  underlying  physical 
mechanisms  will  more  likely  be  made  if  the  problem  is  broken  down  into  various  sub-aspects 
rather  than  if  all  issues  are  addressed  simultaneously.  We  decided  to  advance  simulation 
methodologies  and  understanding  for  related  simplified  model  problems  that  exclude  certain 
aspects  of  the  overall  scenario.  In  particular,  using  both  water  tunnel  experiments  and  direct 
numerical  simulations  (DNS),  we  investigated  3-D  separation  bubbles  on  a  flat  plate.  Separation 
was  induced  by  a  displacement  body  that  was  located  in  close  proximity  to  the  plate. 

The  flow  around  naval  vessels  is  also  characterized  by  high  Reynolds  numbers.  Simulations 
of  such  flows  are  challenging  for  many  reasons:  The  spectrum  of  turbulent  length  scales  spans 
several  orders  of  magnitude.  Even  with  the  most  powerful  supercomputers  it  is  impossible  to 
resolve  all  turbulence  length  scales.  This  motivates  the  development  of  turbulence  models  that 
capture  as  much  resolved  motion  as  permitted  by  the  grid  resolution  and  model  all  remaining 
smaller  scale  turbulent  motion.  For  a  variety  of  engineering  problems  that  involve  large-scale 
unsteadiness,  Reynolds-averaged  Navier-Stokes  (RANS)  does  not  yield  accurate  results  and 
large  eddy  simulation  (LES)  is  often  not  computationally  feasible.  In  particular,  LES  suffers 
from  a  very  restrictive  grid  resolution  requirement  near  walls.  An  idea  pursued  by  many 
researchers  is  to  switch  or  gradually  blend  to  RANS  near  walls.  Speziale  was  among  the  first  to 
propose  a  hybrid  turbulence  modeling  approach  that  combines  the  advantages  of  RANS,  LES, 
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and  DNS  [Speziale  1997],  Ideally,  hybrid  models  blend  DNS,  LES,  and  RANS  such  that  optimal 
use  is  made  of  the  available  computational  assets.  This  requires  that  the  filter  width  is  related  to 
the  spatial  and  temporal  resolution  and  the  employed  turbulence  models  are  state-of-the-art.  In 
Speziale’s  approach  the  “feedback”  from  the  turbulence  model  equations  to  the  Navier-Stokes 
equations  is  scaled  by  a  contribution  function  which  is  dependent  on  the  local  ratio  of  grid  line 
spacing  and  Kolmogorov  length  scale.  This  approach  was  later  called  flow  simulation 
methodology  (FSM)  [Fasel  et  al.  2002,  Israel  2005,  Zhang  et  al.  2000],  Another  example  of  a 
hybrid  model  with  a  fixed  filter  width  is  the  filter-based  RANS  (FBR)  approach  by  Johansen  et 
al.  [2004],  This  model  was  designed  such  that  a  one-equation  sub-grid  stress  (SGS)  model  is 
recovered  in  the  fine  grid  limit.  Because  the  fixed  filter  width  has  to  be  chosen  such  that  it  is 
always  larger  than  the  local  grid  spacing,  the  spatial  and  temporal  resolution  will  not  always  be 
optimal.  This  motivated  Hajjawi  et  al.  [2008]  to  employ  the  FBR  with  a  local  filter  width.  A 
simulation  of  a  turbulent  channel  flow  using  a  hybrid  turbulence  model,  where  a  Smagorinsky 
model  was  switched  to  a  RANS  model  was  carried  out  by  Hamba  [2003].  A  mismatch  in  the 
velocity  profiles  was  observed  at  the  interface  region  between  the  two  approaches.  Turbulent 
channel  flow  simulations  by  Piomelli  et  al.  [2003]  address  this  difficulty.  The  discontinuity  at  the 
interface  is  attributed  to  a  mismatch  of  scales  between  RANS  and  LES.  While  the  turbulence 
model  supports  most  of  the  Reynolds  shear  stress  in  the  RANS  region,  the  resolved  eddies  must 
supply  the  dominant  contribution  in  the  LES  region.  In  the  transition  region  the  eddy-viscosity 
contribution  to  the  mean  shear  is  too  low  while  the  energy-carrying  eddies  have  not  yet  been 
generated.  The  problem  can  be  partially  remedied  through  the  introduction  of  a  “backscatter” 
model  based  on  stochastic  forcing.  For  example.  Batten  et  al.  [2004]  introduced  random  velocity 
fluctuations  obtained  from  a  synthetic  turbulence  model  for  seeding  turbulence  in  their  hybrid 
turbulence  model  simulations  of  a  turbulent  channel  flow. 

Turbulent  channel  and  diffuser  flows  are  good  test  cases  for  such  models.  In  particular,  low 
Reynolds  number  duct  flows  are  frequently  chosen  for  DNS,  RANS,  and  hybrid  model 
validation  efforts.  The  turbulent  asymmetric  diffuser  flow  experiment  at  Stanford  University  is  a 
particularly  well  documented  experiment  [Cherry  et  al.  2008].  In  this  experiment  the  approach 
channel  flow  transitions  fully  before  reaching  the  diffuser  inlet.  The  approach  flow  Reynolds 
number  was  15,380  based  on  hydraulic  diameter  and  bulk  velocity  and  10,000  based  on  channel 
height  and  bulk  velocity.  Earlier  2-D  diffuser  experiments  suffered  from  the  fact  that  an  infinite 
span  cannot  be  realized  in  the  laboratory.  This  issue  is  avoided  for  the  3-D  diffuser  because  of 
the  side  walls.  Using  various  different  hybrid  turbulence  models  we  simulated  the  flow  through 
the  so-called  baseline  diffuser  geometry  for  which  full  flow  field  experimental  data  is  available 
[Cherry  et  al.  2008].  As  these  simulations  are  very  compute  time  intensive  we  decided  to  modify 
and  validate  candidate  hybrid  RANS/LES  models  for  a  turbulent  square-duct  flow  at  a  bulk 
Reynolds  number  of  1 0,000,  because  for  this  Reynolds  number  reference  data  for  validation  are 
available  in  the  literature.  Models  that  fail  to  predict  the  channel  flow  will  likely  also  fail  to 
predict  the  turbulent  separation  in  the  diffuser. 

Turbulent  duct  flows  such  as  the  square-duct  flow  occur  in  many  technical  applications  such 
as  air  ducts  in  buildings  or  ducts  in  aerospace  and  marine  applications.  In  square-duct  flows,  an 
imbalance  of  the  Reynolds-stresses  near  the  comers  results  in  a  secondary  flow  which  transports 
fluid  into  the  comers  [Nikuradse  1926].  With  linear  RANS  models  the  Reynolds-stress 
imbalance  cannot  be  captured  and  the  secondary  flow  does  not  develop.  Predictions  of  the 
secondary  flow  require  non-linear  Reynolds  stress  or  full  Reynolds-stress  models.  The  low- 
Reynolds  number  square-duct  flow  is  accessible  by  DNS.  Huser  and  Biringen  [1993]  carried  out 
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a  DNS  of  a  square-duct  flow  for  a  Reynolds  number  based  on  bulk  velocity  and  hydraulic 
diameter  of  Re  =  1 0,320  (which  corresponds  to  a  Reynolds  number  based  on  mean  friction 
velocity  and  hydraulic  diameter  of  Rex=  600).  The  computational  grid  had  96x101x101  points 
(streamwise  x  cross-flow  directions)  and  dimensions  6.4hxhxh  (length  x  width  x  height),  where 
h  is  the  channel  height.  The  convective  terms  were  discretized  with  fifth-order-accurate  upwind 
biased  finite-differences  and  the  viscous  terms  were  discretized  with  fourth-order-accurate 
central  finite-differences.  The  minimum  correlation  in  the  streamwise  direction  was  at  Ax  =  3.2h, 
indicating  that  the  streamwise  domain  extent  of  Lx  =  6.4h  was  sufficient.  The  comparison  with  a 
RANS  calculation  based  on  Speziale’s  non-linear  k-e  model  revealed  some  significant 
differences  between  the  model  behavior  and  the  DNS  results,  such  as  for  example,  an  under¬ 
prediction  of  the  intensity  of  the  secondary  flow  [Huser  at  al.  1994],  Gavrilakis  [1992]  simulated 
a  square-duct  flow  at  ReT  =  300.  The  computational  domain  for  this  DNS  had  dimensions 
57ihxhxh  with  1000x127x127  points.  The  code  was  second-order-accurate  in  both,  space  and 
time.  In  LES  sub-grid  turbulence  is  modeled  thus  allowing  for  a  lower  grid  resolution  compared 
to  DNS.  Madabhushi  and  Vanka  [1991]  carried  out  a  LES  with  the  Smagorinsky  SGS  model  of  a 
square-duct  flow  at  ReT  =  360  (Re  =  5,810).  The  grid  resolution  was  65x65  in  the  cross-flow 
plane  with  32  Fourier  modes  in  the  streamwise  direction.  The  dimensions  of  the  computational 
domain  were  2nhxhxh.  The  discretization  was  second-order-accurate.  We  employed  DNS  and 
RANS  as  well  as  hybrid  RANS/LES  for  simulating  the  square-duct  flow  with  different  grid 
resolutions.  For  the  hybrid  simulations,  we  utilized  both  the  FSM  and  FBR  based  on  different 
versions  of  the  k-oo  model  by  Wilcox  [2006]. 

This  report  is  organized  as  follows:  In  Section  2  the  scope  of  investigations  is  introduced. 
Section  3  provides  details  on  the  simulation  strategies  and  water  tunnel  experiments.  Results  are 
provided  in  Sections  4-7.  Finally,  in  Section  8  the  results  are  summarized  and  conclusions  are 
drawn. 
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2.  SCOPE  OF  INVESTIGATIONS 


The  main  objective  of  our  numerical  investigations  and  water  tunnel  experiments  was  to  gain 
insight  into  the  physical  mechanisms  governing  3-D  separation.  Of  particular  interest  were  the 
development,  dynamics,  and  interaction  of  energetic  vortical  structures  (spanwise/streamwise) 
that  develop  from  the  different  hydrodynamic  instabilities  that  are  present  in  3-D  separated 
flows. 

We  also  investigated  the  flow  through  the  Stanford  asymmetric  diffuser  by  J.  Eaton  and  co¬ 
workers  [Cherry  et  al.  2008],  which  is  a  challenging  benchmark  for  evaluating  the  performance 
of  different  turbulence  modeling  approaches  such  as  3-D  RANS  and  hybrid  RANS/LES.  Of 
particular  interest  here  was  the  development  of  improved  turbulence  models  for  complex,  non¬ 
equilibrium  turbulent  shear  flows.  As  a  reference,  and  for  validation  of  our  simulation  codes,  we 
also  considered  the  turbulent  square-duct  flow.  This  flow  is  less  challenging  but  better 
documented  than  the  diffuser  flow  and  thus  served  as  an  additional  test  case  for  the  various 
turbulence  models. 

In  summary,  our  approach  consisted  of: 

1)  High-resolution  DNS  of  3-D  separation  bubbles  on  a  flat  plate  (§4) 

2)  Water  tunnel  experiments  of  3-D  separation  bubbles  on  a  flat  plate  (§5) 

3)  Turbulent  square-duct  flow  simulations  (§6) 

4)  Simulations  of  the  turbulent  flow  through  the  Stanford  baseline  asymmetric  diffuser  (§7) 
We  considered  three  different  geometries: 

1)  We  carried  out  simulations  and  experiments  where  we  investigated  3-D  separation 

bubbles  on  a  flat  plate.  Separation  was  induced  through  the  close  vicinity  of  an  asymmetric 

displacement  body. 

2)  For  turbulence  model  validation  we  simulated  the  turbulent  flow  through  a  square-duct. 

3)  We  simulated  the  turbulent  flow  through  the  Stanford  asymmetric  diffuser  which  was 

designed  and  tested  by  J.  Eaton  and  co-workers  at  Stanford  University. 
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2.1  Geometry  1:  Axisymmetric  Displacement  Body  over  Flat  Plate 

In  the  water  tunnel  experiments,  an  axisymmetric  cone-shaped  displacement  body  of 
diameter  D  with  a  spherical  front  end  was  placed  in  close  proximity  above  a  flat  plate  to  generate 
a  separation  bubble  on  the  flat  plate  (see  Fig.  2.1).  Suction  was  applied  through  holes  in  the 
displacement  body  in  order  to  prevent  separation  at  the  aft  section  of  the  displacement  body. 

All  measurements  and  computations  are  for  a  body  diameter  of  D=0. 1  m  and  a  half-opening  angle 
of  20°.  By  adjusting  the  distance  H  of  the  displacement  body  from  the  flat  plate,  the  non- 
dimensional  pressure  distribution  can  be  altered.  The  boundary  layer  thickness  at  separation  can 
be  varied  by  either  moving  the  body  in  the  streamwise  direction  or  by  changing  the  approach 
flow  velocity  Uoo.  For  both  simulations  and  experiments  the  distance  to  the  leading  edge  of  the 
flat  plate  was  kept  constant  at  s  =  0.25m. 


U, 


X 


c> 


►1  axisymmetric  3-D  displacement  body 


H 


flat  plate  3-D  separation  bubble 

Fig  2.1  Axisymmetric  displacement  body  over  flat  plate. 


2.2  Geometry  2:  Stanford  Asymmetric  Diffuser 

We  simulated  the  flow  through  the  asymmetric  3-D  diffuser  by  Eaton  and  co-workers  (Fig. 
2.2)  for  which  high  resolution  velocity  volume  data  is  available  [Cherry'  et  al.  2008].  The 
rectangular  duct  approach  flow  has  cross-flow  dimensions  h><3.33h  (height  x  width)  and  a  length 
that  is  sufficient  to  guarantee  fully  turbulent  flow  at  the  diffuser  inlet.  The  approach  flow 
Reynolds  number  is  10,000  based  on  duct  height  and  bulk  velocity  and  15,380  based  on 
hydraulic  diameter  and  bulk  velocity.  The  diffuser  opening  angles  are  atan(4 —  1 )/ 1 5= 1 1.3deg  and 
atan(4-3.33)/15=2.56deg.  The  cross-sectional  area  at  the  diffuser  exit  is  4x4=16.  The 
experimental  data  shows  turbulent  separation  in  the  diffuser. 
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development  channel 


diffuser  inlet 

h=1cm,  b=3.33cm 
A=3.33cm2 


diffuser  outlet 
4x4cm 
A=16cm2 


Side 


Top 


Fig.  2.2:  Stanford  baseline  diffuser. 


2.3  Geometry  3:  Square-Duct 

In  addition,  we  investigated  the  turbulent  flow  through  a  square-duct  at  a  Reynolds  number 
based  on  hydraulic  diameter  and  bulk  velocity  of  10,000.  Because  this  Reynolds  number  is  close 
to  the  Reynolds  number  of  the  diffuser  approach  flow  we  were  confident  that  turbulence 
modeling  strategies  that  were  successful  for  the  square-duct  flow  would  also  likely  be  adequate 
for  diffuser  flow  simulations. 
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3.  COMPUTATIONAL  AND  EXPERIMENTAL  APPROACH 


3.1  Computational  Methods 

For  our  simulations  we  employed  two  different  numerical  methods: 

(1)  A  high-order-accurate  incompressible  DNS  code  (for  geometry  1) 

(2)  A  high-order-accurate  compressible  multi-domain  finite  volume  code  (for  geometries  2  &  3) 

3.1.1  Method  1 :  Higher-Order-Accurate  Incompressible  Finite  Difference  Code 

For  our  Direct  Numerical  Simulations  (DNS)  of  3-D  separation  bubbles  on  a  flat  plate  we 
employed  a  highly  efficient  and  accurate  incompressible  Navier- Stokes  code  originally 
developed  for  flat-plate  boundary  layer  transition  simulations  [Meitz  2000]  and  later  extended  to 
generalized  curvilinear  orthogonal  coordinates  [Postl  2006]. 

3.1.1  1  Governing  Equations 

In  this  code,  the  3-D  incompressible  Navier-Stokes  equations  are  solved  in  vorticity-velocity 
formulation.  The  equations  are  non-dimensionalized  using  a  reference  velocity  Uref  (e.g.  the 
freestream  velocity)  and  a  reference  length  Lref.  By  taking  the  curl  of  the  momentum  equations, 
which  eliminates  the  pressure  terms,  the  transport  equation  for  the  vorticity  (in  vector  form)  is 
obtained, 

—  +  Vx(«xf)=- V2®  +  VxF,  (3-') 


with  the  vorticity  defined  as 


3  =  -VxV.  (3.2) 

A  volume  force  F.  can  be  added  to  the  right-hand-side  to  introduce  disturbances  into  the  flow. 
From  the  vorticity  <y.  the  velocity  field  V  is  obtained  by  solving  a  vector  Poisson  equation, 

V2V  -V  xco.  (3.3) 

The  velocity  vector  V  has  components  in  the  streamwise,  wall-normal,  and  spanwise  direction 
{u,  v.  and  w),  and  the  components  of  the  vorticity  vector  <y  are  defined  as 

dv  dw  dw  du  and  du  dv  (3  4) 

dz  8y  8x  dz  ‘  8y  dx 


3. 1. 1.2  Discretization 

For  the  time-integration  of  the  vorticity-transport  equations,  a  four-stage  explicit 
Runge-Kutta  scheme  with  fourth-order  accuracy  is  employed.  For  the  spatial  discretization, 
fourth-order  accurate  compact  finite  differences  are  used  in  the  streamwise  and  in  the  wall- 
normal  direction.  Assuming  periodicity  in  the  spanwise  direction,  a  pseudo-spectral  approach  in 
z  is  used  with  a  Fourier  decomposition  of  the  velocity  and  vorticity  into  finite  sine  and  cosine 
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series  with  k  =  [0  ...  K]  Fourier  components.  For  the  flow  variables  u,  v,  and  ox  the  Fourier  series 
has  the  form 

d(x,  y,  z,  t)  =  (x,  y,  /)  +  Z  °k  (*.  y>  0  cos  O'*  z)  ~^Dk(x,y,  t)  sin(y*z), 


with  d=  u,v,  ox  and  Dk  =  t/,  F*,  /2/.  For  the  flow  variables  w,  and  ODy  the  Fourier  series  has 
the  form 

“  3*(x,y,/)cos(r*z),  (3’6) 

A:=l  k=- 1 


d(x,  y,  z,t)=Dk-° (x,  *  /) + Z  *>*  (x,  y,  /)sin(y* z)  +  £  Dk  (x,  y,  t) cos (y k z). 


with  d=w,  a)*,  o)y  and  Dk  =  W*,  Qk,  C2k.  The  spanwise  wave  number,  yk,  is  defined  as 

yk  =  2nkAzk=I,  (3.7) 

where  Xf  1  is  the  spanwise  wavelength  of  the  lowest  spanwise  Fourier  mode,  which  effectively 
is  the  domain  width,  L:,  of  the  computational  domain.  In  the  Fourier  sine  and  cosine  expansions 
Eqs.  (3.5)  and  (3.6),  the  two-dimensional  flow  component  is  represented  by  mode  k  =  0.  The 
symmetric  part  of  the  three-dimensional  flow  field  is  represented  by  modes  1  <  k  <  K,  the 
asymmetric  part  by  modes  —K  <  k  <  -1 .  In  all  our  present  simulations  we  assumed  spanwise 
symmetry,  which  halves  the  number  of  spanwise  modes  required  but  introduces  slip  walls  at  z  = 
L:  and  L-/2  where  w,  ox,  (Oy  =  0.  This  confines  the  movement  of  streamwise  structures  in  the  flow 
by  preventing  them  from  meandering  across  the  symmetry  plane. 

A  key  component  of  the  present  code  is  a  fast  and  parallelized  multi-grid  solver  for 
computing  the  velocity-Poisson  equations.  At  every  computational  step,  disturbances  at  grid 
level  in  the  vorticity  components  are  filtered  out  using  a  fourth-order  compact  filter.  This  allows 
for  simulations  on  grids  that  cannot  resolve  turbulent  structures  all  the  way  down  to  the 
Kolmogorov  length  scale. 


3  1. 1.3  Computational  Setup 

At  the  inflow  boundary,  steady  velocity  and  vorticity  profiles  are  imposed.  To  prevent 
reflections  from  the  outflow  boundary  at  xoui,  the  flow  is  relaminarized  in  a  buffer  domain 
starting  at  Xbuff  using  an  approach  similar  to  that  proposed  by  Kloker  et  al.  [Kloker  1993].  At  the 
free  stream  boundary,  ymax,  the  flow  is  assumed  to  be  irrotational  and  all  vorticity  components 
and  their  derivatives  are  set  to  zero.  Also,  a  Dirichlet  boundary  condition  is  applied  for  the  wall- 
normal  velocity.  The  wall-normal  velocity  was  obtained  from  precursor  calculations  that  are 
described  in  Section  4.1.  At  the  wall,  y  =  0,  no-slip  and  no-penetration  conditions  are  imposed. 
Numerous  simulations  with  different  domain  sizes,  grid  resolutions,  baseflow  and  forcing 
parameters,  etc.,  have  been  performed  over  the  course  of  the  present  research  project.  Details  on 
the  computational  parameters  for  the  different  simulations  are  provided  in  Section  4. 
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3. 1.1. 4  Disturbance  Generation 

Disturbances  are  introduced  through  a  blowing  and  suction  slot  on  the  flat  plate.  The  slot  is 
modeled  by  prescribing  a  wall-normal  velocity  distribution 

vL*)  =  4*)sin(*  )3  (3.8) 

on  the  wall  for  each  spanwise  mode  k.  The  time-dependent  amplitudes  ^**  of  the  spanwise 
modes  are  defined  as 


-O')  =  sin(2^*>/  +  O •  (3.|  2) 

with  frequency,/*,  and  phase  angle,  0k.  The  freestream  effect  of  the  slot  forcing  was  modeled 
assuming  potential  flow  over  an  infinitely  long  flat  plate  with  a  velocity  distribution  at  the  wall 
that  corresponds  to  the  blowing  and  suction  slot.  Using  a  Green’s  function,  the  potential  is  given 
by 


®(.x,y)=-^-  jv/£)log[£-x)2 +/]£/£,  (3  9) 

*,1 

where  xsi  and  xS2  are  the  downstream  and  upstream  boundaries  of  the  blowing  and  suction  slot, 
respectively.  The  change  in  the  freestream  velocity  and  its  derivative  then  become 


v(x,y)  =  —  k(£)— 
In  J  (r- 


-y 


(£■ 


■x)2+y2 


—  tx  y)  =  —  fv 

A'  ,y  J  w^,uy~ 


and 


dy 


2  n 


2tg-*y-y^. 

[({-x)2  +  y2]2 


(3.10) 


With  this  method,  symmetric  or  asymmetric  disturbances  can  be  introduced  at  any  downstream 
location. 


3.1.2  Method  2:  Higher-Order- Accurate  Compressible  Multi-Domain  Finite  Volume  Code 

3. 1.2.1  Governing  Equations 

The  Favre-averaged  compressible  Navier-Stokes  equations  are  solved  in  conservative  form 
on  structured  grids.  A  variety  of  turbulence  models  was  implemented  for  Reynolds-Averaged 
Navier-Stokes  (RANS)  calculations  of  turbulent  flows.  For  the  current  investigations  [Gross  and 
Fasel  2008b,  2009a.  2010],  the  1988,  1998,  and  2006  versions  of  the  k—t o  model  by  Wilcox 
[Wilcox  2006]  were  employed.  The  1998  and  2006  versions  allow  for  improved  predictions  of 
shear  layers.  The  Menter  shear-stress  transport  (SST)  model  [Menter  1994]  which  removes  the 
k-co  models  sensitivity  to  the  free-stream  value  of  the  turbulence  specific  dissipation  [Menter 
1992],  the  k-s  model  [Jones,  1972],  the  Lam-Bremhorst  low-Reynolds  number  k-z  model  [Lam 
1981],  and  the  Spalart-Allmaras  one-equation  turbulence  model  [Spalart  1992]  were  considered 
as  well.  In  the  original  formulations  of  the  various  turbulence  models  addressed  above,  the 
Boussinesq-approximation  (in  the  following  referred  to  as  B-A)  is  employed  for  computing  the 
Reynolds  stresses, 


*V  =—zPkSv+2ftT\S¥- 


-  V 


k  ,k  ij 


(3.11) 


1 1 


where  k  is  the  turbulence  kinetic  energy,  Sv=l/2(Uij+Uj  j)  and  W,j=\ /2(u,j-Ujj)  are  the  strain  and 
shear  rate  tensor,  and  5,y  is  the  Kronecker  symbol.  Repeated  indices  indicate  summations  while 
commas  in  the  subscript  denote  partial  derivatives.  The  B-A  assumes  a  linear  relation  between 
the  strain  tensor  and  the  Reynolds  stress  tensor.  The  B-A  can  therefore  be  categorized  as  a  linear 
eddy  viscosity  model  (LEVM).  Alternatively,  a  second-moment  closure,  the  explicit  algebraic 
stress  model  (EASM)  in  the  form  proposed  by  Rumsey  and  Gatski  [Rumsey  2001], 


*■/  =-^PkS,J  +2Mt 


s,j  -~vk.ksij 


+  a2a4(slkWkJ  +SjjrM)-2a3a4( SlkS,  -±SHSUSV 


(3.12) 


can  be  employed. 

A  large  number  of  different  hybrid  RANS/LES  models  were  proposed  for  closing  the  gap 
between  RANS  and  LES.  We  considered  the  flow  simulation  methodology  (FSM)  [Speziale 
1997,  Fasel  et  al.  2002]  and  the  filter-based  RANS  (FBR)  approach  [Johansen  et  al.  2004].  For 
the  flow  simulation  methodology  (FSM)  the  turbulence  model  contribution  is  scaled  by  a 
“contribution  function”,  f,  which  in  the  original  formulation  by  Speziale  is  a  function  of  the  local 
grid  spacing,  A=max(A*,A>’,Az),  and  the  Kolmogorov  length  scale  (the  characteristic  length  scale 
of  the  dissipating  eddies), 


Lk  = 


VV« 


\£  J 


(3.13) 


Israel  [2005]  suggested  that  a  third  length  scale,  the  “dissipation  length  scale"  (the  typical  length 
scale  of  the  stress-bearing  motion), 


L 


£ 


k 3  2 


£ 


(3.14) 


which  is  larger  than  the  integral  length  scale,  c^kv 2/e,  should  be  considered.  For /=  0  all  scales  of 
motion  are  resolved  and  the  model  contribution  is  zero;  In  the  RANS  limit  for /=  1  the  model 
contribution  is  1 00%  and  all  unsteady  turbulent  motion  is  modeled.  Effectively,  FSM  simulations 
can  be  considered  to  be  locally  DNS,  large  eddy  simulation  (LES),  or  RANS  depending  on  the 
local  and  instantaneous  ratio  of  A/Z,*. 

Two  alternatives  were  suggested:  (1)  Only  the  terms  that  “feed  back”  into  the  Navier-Stokes 
equations  are  scaled.  The  turbulence  model  is  “dragged  along“  and  influenced  only  indirectly 
since,  as  a  result  of  the  scaling,  either  the  mean  flow  is  changed  or  the  flow  becomes  unsteady. 
(2)  All  turbulence  model  terms  including  those  that  feed  back  into  the  turbulence  model 
equations  itself  are  scaled.  For  the  present  results  we  employed  approach  (2). 

Different  variants  of  the  contribution  function  are  available.  For  example.  Zhang  et  al.  [2000] 
proposed  a  modified  Speziale  contribution  function. 


< 

-  exp 

-  0.001  max 

0,— -2 

l  At  )_ 

(3.15) 
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Alternatively,  and  this  is  the  approach  that  we  employed  for  the  results  shown  here,  the 
contribution  function  can  be  computed  from  the  turbulence  energy  spectrum  (Fig.  3.1). 


length  scale  of  stress  Kolmogorov  length  scale  local  grid  resolution 
bearing  motion 

Fig  3.1:  Turbulence  energy  spectra 


Estimates  for  the  unresolved,  ku,  and  total,  k,  turbulence  kinetic  energy  can  be  obtained  by 
integrating  the  turbulence  energy  spectrum,  E(k),  from  Lk  to  LA  and  Lk  to  LE,  e.g. 


2*7Z., 

k=  J  E(K)dK. 

2x11., 

The  ratio  of  ku  and  k  gives  the  contribution  function. 


(3.16) 


(3.17) 


For  Le  <  Lk,  f  is  set  to  1 .  As  E(k)  we  chose  the  von  Karman  [  1 948]  energy  spectrum.  The  length 
scales  were  computed  as  LE  =  f 1  “LEU  and  Lk  =  f1  4LkU. 

For  filter-based  RANS  (FBR)  [Johansen  et  al.  2004]  the  RANS  eddy  viscosity  is  multiplied 
by 


/=  min 


L.s 
1  C  - 

,’i-3  3/ 

k  2 


(3.18) 


where  the  filter  width,  LA,  is  fixed  in  the  original  formulation.  Similar  to  Hajjawi  et  al.  [2008]  we 
decided  to  use  a  local  filter  width  identical  to  the  local  grid  line  spacing.  We  also  set  C3  to  1 . 

For  both  approaches,  FSM  and  FBR.  which  are  described  in  more  detail  in  [Gross  and  Fasel 
2009a, b]  we  attempt  to  capture  “backscatter”.  The  term 
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(3.19) 


pk, 

f 


df  df 

—  +  u  — 

dt  dx 


+  w 


<r 

dz 


is  added  to  the  RHS  of  the  pk-equation  and  subtracted  from  the  RHS  of  the  pe-equation.  In 
addition,  -pu,’(ft+u,fj)  is  added  to  the  RHS  of  the  momentum  equations  for  ft  +  Ujfj  <  0. 
Artificial  velocity  fluctuations,  Uj’,  are  obtained  from  the  synthetic  turbulence  model  proposed  by 
Batten  et  al.  [2004] 

For  the  channel  flow  simulations  the  term, 


1  2.1  4  - 

-puh  f -  =  - pur 

2  b  Dh  Dhy 


(3.20) 


was  added  to  the  RHS  of  the  x-momentum  equation  to  compensate  for  the  streamwise  pressure 
drop.  Here,  Ub  is  the  bulk  velocity  (mass  flux  divided  by  density  and  cross-sectional  area).  The 
hydraulic  diameter,  Dh  =  4A/U,  is  computed  from  the  channel  cross-section.  A,  and 
circumference,  U.  For  turbulent  channel  flow,  the  friction  factor,  f,  can  be  obtained  from 
empirical  relations,  such  as  the  one  given  by  Petukhov, 

/  =  (0.791nRe,-1.64)-2,  '  (3.21) 


where  the  Reynolds  number  is  based  on  bulk  velocity  and  hydraulic  diameter.  The  relationship 
between  bulk  velocity  and  mean  friction  velocity  is. 


/  = 


\U*J 


(3.22) 


We  also  added  volume  forcing  terms  to  the  continuity  and  energy  equations 

p,  + ...  =  a{pre/  -  p)  (3.23) 

{pe\+...  =  apc\Tr'f-T)  (3.24) 

with  a  =  10  4  to  stabilize  density  and  temperature  for  the  channel  flow  simulations  (we 
employed  a  compressible  code). 

3. 1.2. 2  Discretization 

The  governing  equations  are  solved  in  curvilinear  coordinates.  For  robustness  (especially  on 
highly  distorted  grids)  a  finite  volume  method  is  employed.  The  convective  terms  of  the  Navier- 
Stokes  equations  are  discretized  with  a  ninth -order-accurate  upwind  scheme  based  on  a  weighted 
essentially  non-oscillatory  extrapolation  of  the  characteristic  variables  and  the  Roe  scheme 
[Gross  and  Fasel  2002,  2008a].  We  also  considered  the  second-order-accurate  symmetric  total 
variation  diminishing  (TVD)  scheme  by  Yee  [1987].  This  choice  was  made  since  according  to 
Margolin  and  Rider  [2002]  certain  second-order-accurate  upwind  schemes  have  similar 
properties  as  standard  LES  sub-grid  models.  The  convective  terms  of  the  turbulence  model 


14 


equations  are  discretized  with  a  second-order-accurate  discretization  analogous  to  the  symmetric 
TVD  scheme  by  Yee  [1987],  Fourth-  and  second-order-accurate  discretizations,  respectively,  are 
employed  for  computing  the  Navier-Stokes  viscous  terms  and  the  turbulence  equations  diffusion 
terms.  An  implicit  second-order-accurate  Adams-Moulton  method  is  used  for  advancing  the 
governing  equations  in  time.  The  resulting  system  of  equations  is  solved  iteratively  by  a  Newton 
iteration  based  on  a  line  Gauss-Seidel  algorithm.  The  convergence  of  the  implicit  method  is 
monitored  by  considering  the  root  mean  square  (RMS)  of  the  residuals  of  all  cells  for  each 
individual  equation.  For  the  time  dependent  results  shown  in  this  paper  the  equations  were 
advanced  to  the  next  time  step  when  all  RMS  residuals  dropped  below  0.05.  The  code  was 
parallelized  using  the  message  passing  interface  (MPI)  library. 

3. 1.2. 3  Inflow  and  OutflowCconditions 

We  considered  walls  to  be  adiabatic.  The  turbulence  kinetic  energy,  k,  was  set  to  zero  at 
walls.  For  the  k-co  models  the  smooth  wall  boundary  condition  [Wilcox  2006], 


1  Nu 
(o  =  — j — — 
Re'  pAy 


(3.25) 


was  applied  for  co,  where  Ay  is  the  distance  from  the  wall  to  the  center  of  the  first  cell  away  from 
the  wall.  The  parameter  N  was  set  to  1600.  This  value  guarantees  that  for  near  wall  grid 
resolutions  (in  wall  units)  of  y+  =1  the  effective  surface  roughness  k  is  less  than  5,  which 
corresponds  to  a  hydraulically  smooth  surface.  For  the  channel  flow  simulations  periodicity 
conditions  were  applied  in  the  streamwise  direction.  For  the  diffuser  simulations,  at  the  diffuser 
inflow  the  pressure  was  extrapolated  and  all  other  primitive  variables  were  taken  from  a  channel 
flow  simulation.  At  the  diffuser  outflow'  a  non-reflecting  boundary  condition  [Gross  and  Fasel 
2007]  was  applied. 

3. 1.2.4  Non-Dimensionalization 

All  length  scales  were  non-dimensionalized  with  the  channel  height,  h,  which  for  the  square- 
duct  is  identical  to  the  hydraulic  diameter,  I\  -  4A/U  =  4hV4h  =  h.  Velocities  were  non- 
dimensionalized  with  the  bulk  velocity.  Ub,  based  on  the  nominal  given  bulk  Reynolds  number, 

hydraulic  diameter,  and  kinematic  viscosity.  Alternatively,  the  average  friction  velocity,  «r ,  can 

be  employed  for  non-dimensionalization.  The  desired  bulk  velocity  in  our  square-duct  flow 
simulations  was  Reb  =  10,000.  From  Petukhov’s  relationship  a  friction  factor  of  f  =  0.031479  is 
obtained.  The  average  friction  velocity  then  becomes 


Mr  = 


0.06273wa 


(3.26) 


and  the  Reynolds  number  based  on  average  friction  velocity  and  hydraulic  diameter  is  ReT  = 
627.3.  With  the  average  friction  velocity, 


+ 

y 


=  Rer  =  627.3— ,  and 
K  Dh  Dh 


(3.27) 
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+  _  u  _  Re*  u 
Ur  Rer  ub 


u 

=  15.94  — 

ub 


(3.28) 


are  obtained.  Viscosities  were  normalized  with  the  reference  dynamic  viscosity.  Turbulence 
kinetic  energy  and  Reynolds-stresses  were  non-dimensionalized  with  the  bulk  velocity  squared. 

3.1.3  Post-Processing  Tools 

3. 1.3. 1  Vortex  Identification 

For  visualizing  vortical  structures  we  employed  the  Q-criterion  by  Hunt  [1988], 


Q=  mWijWij-MiStjStj , 


(3.29) 


which  identifies  vortical  structures  by  representing  regions  inside  the  flow  where  rotation,  WtJ 
dominates  strain,  Sy. 

3. 1.3. 2  Auto-Correlation 

Discrete  time-dependent  velocity  data,  un,  can  be  auto-correlated. 


(3.31) 


for  m  =  1 ,. .  ,,N  where  u'  =  u  -  u  is  a  disturbance  velocity. 


(3.32) 


is  the  variance,  and  the  index  n  +  m  -  1  is  taken  modulus  N. 


3  1.3.3  Energy  Spectrum 
Energy  spectra. 


(3.33) 


dm 


can  be  computed  from  velocity  data,  un,  where 


(3.34) 
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(3.35) 


2  N 

h m  sinfo.0- 

For  a  linear  wave  number  (k  =  2n/At)  distribution 

m- 1  /  x 

=  *i  +TT7^ 

M  - 1 

^2-^1 
dm  M  - 1 

is  obtained. 


(3.36) 

(3.37) 
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3.2  Experimental  Facility 
3.2.1  Water  Tunnel 


The  Hydrodynamics  Laboratory  of  the  Department  of  Aerospace  and  Mechanical 
Engineering  at  the  University  of  Arizona  is  a  state-of-the-art  research  facility  that  houses  three 
water  tunnels. 


HOAtycomb  ortv« 


Fig.  3.2:  Sketch  of  closed  water  tunnel  at  the  Hydrodynamics  Laboratory  at  the  Aerospace  and 

Mechanical  Engineering  Department. 

Our  experiments  were  conducted  in  the  “high-speed”  closed  surface  water  tunnel.  This  tunnel 
has  a  test  section  with  dimensions  0.5  x  0.7m.  The  maximum  flow  velocity  in  the  test  section  is 
1.34  m/s.  A  side  view  of  the  tunnel  is  given  in  Fig.  3.2.  A  7.5  kW  electric  motor  drives  a  four- 
bladed  turbine.  The  pitch  of  the  turbine  blades  can  be  adjusted  for  altering  the  velocity  range. 
The  water  is  fed  back  through  a  pipe  in  the  basement.  A  0.2m  wide  honey  comb  section  is  used 
for  conditioning  the  flow  (flow  straightening,  turbulence  management).  Some  key  parameters  are 
listed  in  Tab.  3. 1 .  A  chart  showing  the  velocity  range  for  two  different  turbine  blade  pitch  angles 
is  provided  in  Fig.  3.3. 


Tost  Section  Size 

4.5  m  x  0.5  m  x  0.7  m 

Velocity  Range  slow 

0.03  m/s  -  0.58  in/s 

Velocity  Range  high 

0.08  m/s  -  1.34  in/s 

Water  Volume 

56.000  liters 

Tab.  3.1:  Specifications  of  closed  water  tunnel. 
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Tunnel  Speed 

•  WT  setting  25  degr  turbine  blade  angle 
®  WT  setting  1  degr  turbine  blade  angle 


Re_d 
d=0 . lm 


Fig.  3.3:  Tunnel  speed  over  turbine  revolutions  per  second  for  different  turbine  blade  pitch  angles. 
Velocity  data  was  obtained  from  PIV  measurements. 


322  Flow  Diagnostics 

In  general,  dye  flow  visualization  is  less  involved  and  was  therefore  employed  for  fast  and 
qualitative  flow  diagnostics,  which  provided  considerable  insight  into  the  fluid  dynamics.  We 
utilized  both,  red  and  blue  dye,  which  was  delivered  into  the  flow  through  holes,  and 
phosphorescent  dye,  which  was  injected  with  a  needle.  The  advantage  of  phosphorescent  dye  is 
that  it  allows  for  flow  visualizations  in  a  plane  that  is  illuminated  by  a  laser. 


interrogation 

window  ,ma9f 

plane 


Fig.  3.4:  Sketch  of  PIV  system. 
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We  also  employed  Particle  Image  Velocimetry  (PIV)  for  quantitative  flow  measurements.  In 
general,  PIV  requires  a  high  power  laser,  an  optical  arrangement  to  widen  the  laser  beam  into  a 
light  sheet  (cylindrical  lens),  and  a  camera  (Fig.  3.4).  Velocity  fields  are  extracted  from  the 
movement  of  particles  seeded  into  the  flow.  The  particles  are  almost  neutrally  buoyant,  and  are 
therefore  generally  assumed  to  be  accurately  convected  by  the  flow.  The  measurement  area 
within  the  flow  field  is  defined  by  the  position  of  a  laser  sheet  and  the  optical  opening  angle  of 
the  cameras.  Using  two  short  duration  laser  pulses,  a  double-exposure  of  the  flow  field  is 
obtained.  The  velocity  vectors  are  extracted  by  performing  a  mathematical  correlation  analysis  of 
the  two  separate  frames.  The  PIV  interrogation  process  is  repeated  until  all  required  velocity 
information  is  extracted  from  the  snapshots. 


Fig.  3.5'  Photo  of  PIV  system  including  cameras,  laser,  PIV  computer  and  traverse  system  on  the  left. 

Seeded  flow  and  laser  beam  with  displacement  body  in  test  section  on  the  right. 

A  PIV  system  manufactured  by  LaVision  was  employed  for  our  investigations  (funded  by  a 
DURIP  grant  from  ONR).  This  system  consists  of  a  double  pulsed  NdrYAG  laser  (class  4),  two 
high  quality  digital  cameras  and  a  Dell  quad-core  computer  for  post-processing  (Fig.  3.5).  The 
Imager  ProX2M  cameras  with  Nikon  lenses  have  a  maximum  frame  rate  of  30Hz  at  full 
resolution.  We  designed  and  built  a  traverse  system  that  uses  computer-controlled  stepper  motors 
for  moving  the  laser  and  cameras  in  downstream  direction.  The  traverse  system  allows  us  to  map 
out  the  entire  volume  which  contains  the  three-dimensional  separation  bubble  in  fine  slices.  For 
the  results  shown  in  this  report  the  laser  beam  was  pointed  at  a  mirror  located  in  the  rear  of  the 
test  section  for  generating  laser  sheets  in  spanwise  planes  with  a  thickness  of  approximately 
2mm.  Two  cameras  are  used  to  detect  the  particle  movements  within  the  plane.  FIollow  glass 
spheres  with  a  diameter  of  approximately  10  microns,  which  are  neutrally  buoyant,  were  used  as 
tracer  particles. 
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4.  SEPARATION  BUBBLE:  DIRECT  NUMERICAL  SIMULATIONS 


The  objective  of  our  direct  numerical  simulations  (DNS)  and  water  tunnel  experiments  was 
to  investigate  in  detail  the  structure  and  dynamics  of  three-dimensional  (3-D)  separation  bubbles. 
In  this  section,  results  from  our  simulations  are  presented.  The  water  tunnel  experiments  are 
discussed  in  section  §5. 


flat  plate 


Fig  4.1(a):  Schematic  of  set-up  with  displacement  body  and  computational  domain 


Fig.  4.1(b):  Schematic  of  computational  domain  for  DNS. 


4.1  Precursor  Calculations 

Our  computational  approach  involves  two  steps: 

1 .  A  precursor  calculation  of  the  inviscid  flow  around  the  displacement  body. 

2.  A  high-resolution  Direct  Numerical  Simulation  (DNS)  of  the  separated  flow  region  on 
the  flat  plate  with  freestream  boundary  conditions  that  are  extracted  from  the  precursor 
calculation. 

The  purpose  of  the  precursor  calculations  is  to  obtain  freestream  boundary  conditions  for  the 
highly  resolved  DNS  and  to  study  the  effect  of  the  geometric  parameters  on  the  potential 
flowfield.  The  calculations  are  based  on  the  method  described  by  Zedan  and  Dalton  [1978]  for 
computing  the  potential  flow  around  a  body  of  revolution  with  arbitrary  shape.  Corresponding  to 
the  body  shape,  linear  source  and  sink  distributions  on  the  axis  of  rotation  are  determined  and 
superposed  on  a  uniform  flow.  The  source  and  sink  distributions  are  found  by  solving  a  system 
of  linear  equations,  which  is  obtained  from  the  conditions  that  (i)  the  streamfunction  is  zero  at  a 
certain  number  of  control  points  on  the  body  surface,  and  (ii)  the  net  mass  flux  of  the  entire 
source  and  sink  distribution  is  zero.  To  account  for  the  presence  of  the  flat  plate,  a  mirror  image 
is  generated  and  added  to  the  original  solution.  The  change  of  the  original  body  shape  due  to  the 
mirroring  was  found  to  be  negligible.  For  validation  purposes,  the  wall-normal  velocity  profiles 
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obtained  with  this  procedure  were  compared  to  results  of  inviscid  and  viscous  simulations  of  the 
entire  flow  field  (flat  plate  and  displacement  body)  using  the  commercial  Navier-Stokes  solver 
Fluent.  The  profiles  extracted  from  the  potential  flow  field  and  from  the  inviscid  Fluent 
simulations  agree  very  well  and  are  also  in  very  good  agreement  with  the  profiles  extracted  from 
the  viscous  Fluent  simulations  (Fig.  4.2).  Similarly,  the  influence  of  “neighboring”  displacement 
bodies  due  to  the  periodic  boundary  condition  of  the  DNS  was  accounted  for  by  repeating  and 
mirroring  the  solution  in  the  spanwise  direction. 


(ft)  Profiles  of  wall-normal  velocity  in  the  mid-span  (b)  Profiles  of  wall-normal  velocity  in  spanwise  direct  ion 
plane  at  the  downstream  location  of  the  suction  peak. 


Fig.  4.2:  Profiles  of  wall-normal  velocity  at  wall-normal  location  ymax=  0.048m,  extracted  from  potential 
flow  and  Fluent  precursor  calculations  for  different  Reynolds  numbers  Ren  (for  the  baseline  case  with 
H-lm). 


4.2  Simulation  Setup 

The  computational  domain  for  the  high-resolution  DNS  is  shown  in  Fig.  4.1b.  The  inflow 
boundary  is  located  at  a  distance  Xi=0.1m  from  the  leading  edge  of  the  flat  plate.  The  domain 
extends  a  distance  ymax  from  the  surface  of  the  flat  plate.  It  spans  from  the  midspan  plane  at  z  =  0 
to  z  =  z.max.  For  all  simulations  flow  symmetry  with  respect  to  the  body  center  was  assumed.  In 
addition,  the  flow  was  assumed  to  be  periodic  in  the  spanwise  direction.  The  distance  of  the 
displacement  body  from  the  flat  plate  leading  edge  was  s  =  0.25m,  or  0.15m  measured  from  the 
inflow  boundary  of  the  computational  domain.  The  spanwise  variation  of  the  freestream  velocity 
at  the  inflow  boundary  due  to  the  upstream  influence  of  the  displacement  body  was  less  than 
three  percent.  The  streamwise  domain  length  including  the  outflow  was  0.640m  for  the  unsteady 
cases  and  0.768m  for  the  steady  cases.  The  outflow  boundary  was  thus  at  X2  =  0.740m  and 
0.868m,  respectively.  Normalized  with  the  displacement  body  diameter,  the  streamwise  extent  of 
the  computational  domain  measured  from  the  center  of  the  displacement  body  was  1.5D  in 
upstream  direction  and  3.9D  in  downstream  direction. 

The  height  ymax  of  the  computational  domain  was  varied  between  0.048m  and  0.068m.  To 
investigate  the  influence  of  the  spanwise  domain  width  on  the  flow  (Sec.  4.3.1),  domain  widths 
from  0.24m  to  1 .28m,  i.e.  from  W  =  1 .2  to  W  =  12.8  body  diameters,  were  investigated. 
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The  grid  for  the  steady  simulations  (i.e.  for  simulations  where  the  Reynolds  number  is  low 
enough  so  that  the  flow  remained  steady)  had  a  resolution  of  lxlO  'm  in  the  streamwise 
direction.  In  the  wall-normal  direction,  129  or  161  points  were  used  depending  on  ymax  with  a 
mild  exponential  grid  point  clustering  near  the  wall  to  adequately  resolve  the  boundary  layer  on 
the  flat  plate  as  well  as  the  developing  shear  layer  and  vortices.  The  stretching  was  adjusted  such 
that  the  near  wall  resolution  for  the  different  ymax  was  similar  for  the  different  cases  and  about 
half  of  the  grid  points  were  located  in  the  lower  30%  of  the  domain.  The  wall  normal  grid 
resolution  at  the  wall  was  about  1.2xl0_4m.  The  boundary  layer  and  separation  bubble  were 
resolved  with  30—40  points.  A  higher-resolution  grid  was  obtained  by  doubling  the  number  of 
grid  points  in  both  the  streamwise  and  the  wall-normal  direction,  resulting  in  a  streamwise 
resolution  of  0.5x1 0-,m  and  a  near-wall  grid-line  spacing  of  0.6X10  4m.  In  the  spanwise 
direction,  between  61  and  321  modes  were  employed,  which  corresponds  to  193  and  1025 
collocation  points,  respectively,  across  the  full  physical  spanwise  domain  width.  This  results  in  a 
spanwise  resolution  of  1 .25x1 0  3m  for  the  lower-resolution  grid  and  0.625x10  3m  for  the  higher- 
resolution  grid.  For  the  unsteady  simulations,  we  utilized  the  higher-resolution  grid  with  257 
points  in  the  wall-normal,  1281  points  in  the  streamwise  and  up  to  641  modes  in  the  spanwise 
direction  (337  million  points  total).  The  minimum  grid  resolution  in  wall  units  in  the  wake  was 
x+  ~  1 3,  y+  ~  1 .9,  and  z+  ~  9. 

4.3  Results 

To  gain  a  basic  understanding  of  the  topology  of  the  3-D  separation,  several  steady  flow 
cases  were  investigated  first.  In  this  chapter,  the  flow  topologies  for  these  steady  flow  cases  as 
well  as  for  the  time-averaged  flowfields  obtained  from  simulations  at  larger  Reynolds  numbers, 
where  the  flow  is  unsteady,  are  presented.  The  unsteady  features  of  the  latter  are  described  in 
Sec.  4.4  and  the  flow  structures  for  all  cases  are  discussed  in  more  detail  in  Sec  4.5. 

The  flow  topologies  were  analyzed  based  on  the  limiting  streamline  pattern  on  the  flat  plate 
and  streamlines  in  the  symmetry  plane  (Figs.  4.4-4.8).  The  limiting  streamlines  were  overlayed 
with  contours  of  spanwise  vorticity  on  the  plate  and  streamwise  velocity  in  the  symmetry  plane 
(z  =  0)  to  indicate  the  reverse  flow  region,  which  corresponds  to  negative  values  of  vorticity  and 
streamwise  velocity,  respectively.  As  suggested  by  Pauley  (1994),  the  limiting  streamline 
patterns  were  computed  from  the  wall-tangential  velocity  components  at  the  wall  nearest  grid 
points.  In  general,  this  worked  very  well  and  the  constructed  streamline  patterns  are  in  very  good 
agreement  with  the  flowfield  as  indicated  by  the  wall  vorticities. 

4.3.1  Effect  of  Domain  Width 


As  discussed  in  Sec.  3.1.1,  in  the  DNS  code  periodicity  in  the  spanwise  direction  is  assumed. 
Therefore,  the  first  step  was  to  analyze  the  effect  of  the  domain  width  (W)  on  the  potential  flow 
field  (Fig.  4.3)  and  on  flow  separation  (Fig.  4.4).  The  spanwise  periodicity  implies  that  the 
simulated  scenario  is  not  a  single  separation  bubble,  but  an  infinite  row  of  separation  bubbles,  or, 
in  other  words,  an  infinite  row  of  displacement  bodies.  The  test  section  side  walls  of  the  water 
tunnel  can  also  be  regarded  as  symmetry  planes,  which  for  potential  flow  can  be  modeled  by 
adding  images  of  the  displacement  body  in  the  spanwise  direction.  The  spanwise  periodicity  in 
the  DNS  is  therefore  consistent  with  the  accompanying  experiments.  A  domain  width  of  w  = 
0.64m  or  W  =  6.4  (in  displacement  body  diameters)  was  chosen  as  the  baseline  case,  as  this  is 
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close  to  the  spanwise  width  of  the  test  section  in  the  water  tunnel  (Sec.  3.2)  and  because  this 
results  in  an  inviscid  pressure  distribution  on  the  flat  plate  that  closely  resembles  that  for  a  single 
displacement  body  (Fig.  4.3). 

The  domain  width  mostly  affects  the  spanwise  Cp  gradient  with  little  effect  on  the  stream  wise 
pressure  gradient  (Tab.  4.1).  and,  therefore,  allows  us  to  investigate  the  influence  of  the  spanwise 
pressure  distribution  on  the  separation  topology. 


z[d] 


x  [m] 


(a)  Effect  of  distance  to  "neighboring"  displace¬ 
ment  bodies  on  the  spanwise  Cp  distribution. 
Profiles  shown  are  plotted  at  the  downstream  lo¬ 
cation  of  minimum  Cv. 


(b)  Illustration  of  the  shape  of  C'p  for 
reference  case 

d  =  0.lm,H  =  1.0,1V  =G.4 


Fig.  4.3:  Non-dimensional,  inviscid  pressure  distribution  Cp  on  flat  plate. 

The  analysis  of  the  limiting  streamline  patterns  (Fig.  4.4)  reveals  that  the  width  and,  most 
importantly,  the  topology  of  the  separation  are  greatly  affected  by  the  spanwise  extent  of  the 
domain.  For  all  domain  widths  from  W  =  2.4  to  W  =  6.4,  the  topology  of  the  flow  for  Ren  = 
7,500  and  FI  =  1 .0  shows  one  saddle  point  associated  with  flow  separation  at  the  upstream  end  of 
the  separation  bubble.  From  this  saddle  point  the  line  of  primary  separation  leads  to  a  focus  on 
either  side  of  the  symmetry  plane  and  a  saddle  point  of  attachment  at  the  downstream  end  of  the 
separated  region  (Figs.  4.4a, b). 

When  the  domain  width  is  increased,  the  saddle  point  of  attachment  splits  in  two.  These  two 
separate  saddle  points  move  away  from  the  symmetry  plane,  where  a  node  now  replaces  the 
original  single  saddle  point.  Along  with  this,  a  region  develops  adjacent  to  the  symmetry  plane 
inside  and  downstream  of  the  separation  bubble,  where  the  limiting  streamlines  originate  from 
this  node  of  attachment  (Figs.  4.4c, d).  For  the  two  widest  domains  that  we  considered,  the 
topology  of  the  separation  bubble  remains  unchanged.  This  is,  because  the  inviscid  pressure 
distribution  on  the  flat  plate  (Fig.  4.3a)  for  the  maximum  simulated  domain  width  W  =  12.8  is 
essentially  identical  to  the  pressure  distribution  for  a  domain  with  infinite  spanwise  extent.  Very 
similar  results  are  obtained  for  a  spanwise  domain  width  of  W  =  9.6.  Although  the  difference  in 
the  spanwise  inviscid  pressure  distribution  for  W  =  2.4  and  W  =  6.4  is  much  larger  than  the 
difference  seen  for  W  =  6.4  and  W  =  9.6  (Fig.  4.3a),  the  change  in  topology  is  less  significant  for 
W  =  2.4  and  W  =  6.4  than  for  W  =  6.4  and  W  =  9.6.  The  streamlines  in  the  symmetry  plane 
exhibit  the  same  general  pattern  for  all  cases  and  are  not  affected  by  the  domain  width. 
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index 

H 

w 

Rcd 

RCQ  ,.ir  p 

ACpin„/  A  A' 

ACpjni,/AZ 

i 

1.00 

3.2 

7.500 

102 

0.197 

0.118 

ii 

1.00 

4.8 

7.500 

100 

0.198 

0.096 

iii 

1.00 

G.4 

7.500 

99 

0.197 

0.076 

iv 

1.00 

9.6 

7.500 

98 

0.200 

0.053 

V 

1.00 

12.8 

7,500 

98 

0.199 

0.040 

vi 

1.02 

64 

7,500 

103 

0.185 

0.071 

vii 

1.05 

64 

7,500 

107 

0.167 

0.065 

viii 

1.10 

64 

7.500 

117 

0.142 

0.057 

ix 

1.15 

6.4 

7,500 

— 

0.121 

0.049 

X 

1.00 

6.4 

2.500 

0.197 

0.076 

xi 

1.00 

6.4 

5,000 

88 

0.197 

0.076 

xii 

1.00 

6.4 

15,000 

132 

0.197 

0.076 

index 

L  bubble 

^  bubble 

H  bubble 

Topology 

i 

1.49 

0.52 

0.044 

2  SP,  2  F 

ii 

1.53 

0.54 

0.059 

2  SP.  2  F 

iii 

1.58 

0.58 

0.063 

2  SP.  2  F 

iv 

1.60 

0.64 

0.065 

3  SP.  2  F,  1  N 

V 

1.61 

0.68 

0.065 

3  SP,  2  F,  1  N 

vi 

1.26 

0.52 

0.045 

2  SP.  2  N 

vii 

1.14 

0.50 

0.037 

1  SP,  1  N 

viii 

0.62 

0.30 

0.010 

1  SP,  1  N 

ix 

no  separation 

X 

— 

— 

no  separation 

xi 

0.81 

0.38 

0.028 

1  SP,  1  N 

xii 

sjI.82 

0.88 

0  069 

3  SP.  2  F.  1  N 

Tab  4.1:  Summary'  of  simulation  parameters  (top)  and  separation  bubble  characteristics  (bottom.  SP  = 
saddle  point,  F  =  focus,  N  =  node). 


4.3.2  Effect  of  Displacement  Body  Distance  from  Flat  Plate 

The  distance  of  the  displacement  body  from  the  flat  plate  affects  the  pressure  distribution  on 
the  flat  plate.  While  the  spanwise  domain  width  mostly  affects  the  spanwise  pressure 
distribution,  the  parameter  H  (defined  as  the  displacement  body  distance  from  the  flat  plate)  was 
found  to  affect  both  the  streamwise  and  spanwise  pressure  gradient. 

For  the  simulations  shown  in  Fig.  4.5,  the  displacement  body  was  successively  moved  farther 
away  from  the  flat  plate,  thus  weakening  the  pressure  gradients  acting  on  the  boundary  layer  on 
the  flat  plate.  When  comparing  the  separation  topology  for  increased  distances  H  relative  to  the 
reference  case  with  H  =  1 .0  in  Fig.  4.4b  (Reo  =  7,  500,  W  =  6.4),  it  can  be  seen  that  as  H  is 
increased,  the  saddle  point  of  reattachment  remains  in  the  symmetry  plane,  whereas  the  saddle 
point  of  separation  is  replaced  by  a  node  of  separation  and  the  two  foci  vanish.  This  topology 
appears  to  remain  stable  until  the  pressure  gradient  becomes  too  weak  and  the  separation 
disappears  altogether  (Fig.  4.5d). 

A  closer  look  at  the  streamline  patterns  for  H  =  1.02  and  H  =  1.05  (see  Fig.  4.6)  reveals  a 
slight  but  important  difference,  which  is  seen  most  clearly  in  the  color  contours  of  streamwise 
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vorticity  near  the  line  of  separation.  For  H  =  1 .05  the  line  of  zero  streamwise  vorticity  (indicated 
by  the  white  color)  is  found  downstream  of  the  line  of  separation,  whereas  for  H  =  1.02  it  is 
located  upstream  of  the  line  of  separation.  This  difference  is  crucial,  as  it  indicates  that  for  H  = 
1.02  the  streamlines  converging  on  the  line  of  separation  lead  away  from  the  symmetry  plane 
(same  as  for  H  =  1 .0),  while  for  H  =  1 .05  they  point  towards  the  symmetry  plane.  Consequently, 
the  singular  point  in  the  symmetry  plane  from  which  the  line  of  separation  emanates,  must  be  a 
saddle  point  for  H  =  1 .02,  and  a  node  for  H  =  1 .05.  Since  the  point  of  attachment  is  also  clearly  a 
saddle  point  (in  both  cases),  there  must  be  nodes  or  foci  at  the  outer  ends  of  the  line  of  separation 
for  H  =  1 .02,  according  to  the  topological  rules  first  proposed  by  Lighthill  [1963], 

As  nothing  suggests  a  “spiral”  motion  close  to  the  end  of  the  line  of  separation,  the  points  are 
identified  as  a  pair  of  nodes.  In  the  cases  (with  separation)  where  no  foci  are  present  (Figs.  4.5a- 
c),  the  line  of  separation  has  no  definite  end  point.  (This  fact  necessitates  an  alternative  definition 
for  the  width  of  the  separated  flow  region  [Jacobi  et  al.  2008]).  The  lack  of  definite  end  points  of 
the  line  of  separation  also  represents  a  unique  property  of  3-D  separation  called  “open 
separation”,  which  is  discussed  in  more  detail  in  Sec.  4.3.4.  The  general  pattern  of  the 
streamlines  in  the  symmetry  plane  is  again  not  affected  by  the  changes  seen  in  the  limiting 
streamline  pattern  on  the  flat  plate. 


1  Upon  close  examination  of  the  streamlines  for  H  =  1  02,  one  can  see  that  near  the  outer  end  of  the 
line  of  separation  the  streamline  direction  is  not  in  agreement  with  the  flow  direction  suggested  by  the 
wall-vorticity.  The  streamwise  vorticity  indicates  flow  away  from  the  symmetry  plane  (red),  where  the 
streamlines  turn  towards  the  symmetry  plane.  The  line  of  separation  still  ends  in  a  node,  as  indicated  both 
by  the  streamline  behavior  and  the  wall-vorticity.  The  discrepancy  can  be  attributed  to  the  fact  that  the 
streamline  pattern  was'  computed  from  flow  data  at  a  finite  distance  from  the  plate,  which  may  differ  from 
the  flow  field  “right  on"  the  surface.  Such  a  difference  has  been  observed for  this  case  only. 
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Fig  4.4:  Effect  of  spanwise  domain  width ,  W,  on  separation  topology. 

Shown  are  streamline  patterns  and  color  contours  of  streamwise  velocity  in  the  symmetry >  plane 
(z=0,  on  top),  and  limiting  streamlines  and  color  contours  of  spanwise  vorticity  on  the  plate  (y=0, 
below).  Blue  indicates  reverse  flow. 
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F/g  4.5.  Effect  of  distance  of  displacement  body  from  flat  plate ,  H,  on  separation  topolog}\ 

Shown  are  streamline  patterns  and  color  contours  of  streamwise  velocity  in  the  symmetry  plane 
(z=0,  on  top),  and  limiting  streamlines  and  color  contours  of  spanwise  vorticity  on  the  plate  (y=0 , 
below)  Blue  indicates  reverse  flow . 
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Fig  4.6:  Enlarged  views  of  the  separated  region.  Shown  are  limiting  streamline  patterns  overlayed  with 
color  contours  of  streamwise  vorticity  on  the  flat  plate  (y  =  0).  Blue  indicates  flow  going  "  upwards  ”  (in 
positive  z-direction)f  red  indicates  flow  going  “downwards”  (in  negative  z-direction). 


4.3.3  Effect  of  Reynolds  Number 

A  wide  range  of  Reynolds  numbers  Reo  (based  on  displacement  body  diameter  and 
freestream  velocity)  between  2,500  and  30,000,  were  investigated  for  a  fixed  spanwise  domain 
width  W  =  6.4  and  distance  of  the  displacement  body  from  the  flat  plate  H  =  1.0.  The  separated 
flow  is  steady  up  to  a  Reynolds  number  of  10,000,  it  is  unsteady  for  Ren  =  1 1,000,  15,000,  and 
30.000.  For  the  unsteady  cases,  the  flow  topologies  were  deduced  from  the  time-averaged 
flowfields. 

Fig.  4.7  shows  the  streamline  pattern  in  the  symmetry'  plane  and  on  the  flat  plate  for  the 
steady  cases.  For  Ren  =  5,000  (Fig.  4.7b)  the  flow  topology  is  the  same  as  for  some  of  the  cases 
where  the  height  of  the  displacement  body  was  increased  (Figs.  4.5b,c).  The  limiting  streamline 
pattern  on  the  flat  plate  only  contains  one  node  of  separation  and  one  saddle  point  of 
reattachment,  and  the  line  of  primary'  separation  shows  the  same  “open  separation  behavior  ’  as 
described  in  Sec.  4.3.2  for  H  =  1.05.  (For  Ren  =  2,500  and  H  =  1.0,  as  well  as  for  ReD  =  7,500 
and  H  =  1.15,  the  flow  did  not  separate  from  the  flat  plate.) 

The  streamline  patterns  for  increased  displacement  body  distances  H  and  decreased  Reynolds 
numbers  Reo  are  similar.  Varying  the  Reynolds  number  Reo,  however,  does  not  directly  affect 
the  non-dimensional  inviscid  pressure  distribution  Cp  mv,  as  it  is  independent  of  the  freestream 
velocity.  Thus,  for  reproducing  a  particular  experiment,  both  Reynolds  number  and  Cp  have  to  be 
matched.  Nevertheless,  the  simulations  show  that  it  is  possible  to— at  least  qualitatively — 
duplicate  the  flow  behavior  at  a  different  Reynolds  number,  if  the  Cp  distribution  imposed  on  the 
flat  plate  is  adequately  adjusted. 

For  the  unsteady  cases,  the  time-averaged  flowfields  shown  in  Fig.  4.8a-c  indicate  topologies 
that  to  some  extent  differ  from  the  reference  case  with  Reo  =  7,500  and  W  =  6.4,  but  also  share 
some  commonalities  with  the  Reo  =  7,500  cases  with  large  domain  widths:  For  example,  the  line 
of  primary  separation  connects  a  saddle  point  of  separation  with  two  focal  points,  although  the 
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foci  are  much  stronger  for  the  higher  Reynolds  numbers.  Also,  the  width  of  the  bubble,  as 
determined  by  the  spanwise  extrema  of  the  line  of  primary  separation,  is  considerably  larger. 
While  the  saddle  point  of  separation  is  seen  to  move  upstream  with  increasing  Reynolds  number, 
the  location  of  the  foci  remains  almost  unchanged,  compared  to  the  reference  case. 

The  attachment  regions  of  the  time-averaged  flowfields  exhibit  a  topology  that  is  not  easily 
interpreted,  but  qualitatively  similar.  A  node  of  attachment  is  located  in  the  symmetry  plane  with 
a  saddle  point  on  either  side  (Fig.  4.8a-c).  The  node  of  attachment  is  necessitated  by  the  theorem 
governing  the  number  of  singular  points,  provided  one  accepts  the  notion  that  the  singular  points 
away  from  the  symmetry  plane  are  saddle  points  [Lighthill  1963]. 

Compared  to  the  topology  of  a  separation  bubble  at  Ren  =  7,500  (W  =  12.8,  Fig.  4.4d),  which 
displays  an  equal  number  and  arrangement  of  singular  points,  the  node  of  attachment  in  the 
unsteady  cases  appears  to  be  rotated  by  90  degrees,  and  the  streamline  patterns  near  the 
symmetry  plane  and  downstream  of  the  point  of  reattachment  differ.  Figure  4.8  also  shows  that 
at  least  one  line  of  secondary  separation  emanates  on  either  side  of  the  symmetry  plane  from  the 
saddle  points  downstream  of  the  separation  bubble.  Note  that  for  all  the  cases  which  we 
investigated  this  is  the  only  indication  of  a  secondary  separation. 

The  difference  in  the  separation  topology  between  the  unsteady  and  steady  cases,  which 
shows  up  in  the  streamline  patterns  in  the  symmetry  plane,  can  also  be  seen  in  Fig.  4.10. 
Whereas  for  the  steady  reference  case  the  shear  layer  in  the  symmetry  plane  stays  detached  from 
the  wall  (Fig.  4.10a),  for  the  unsteady  cases  the  flow  reattaches  for  the  time-averaged  flow  field 
(Fig.  4.10b-d).  This  is  discussed  in  more  detail  in  Sec.  4.4. 
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F/g  4.  7:  Effect  of  Reynolds  number  Reo  on  separation  topology  (steady  cases). 

Shown  are  streamline  patterns  and  color  contours  of  streamwise  velocity  in  the  symmetry?  plane 
(z=Of  on  top),  and  limiting  streamlines  and  color  contours  of  spanwise  vorticity  on  the  plate  (y=0, 
below).  Blue  indicates  reverse  flow. 
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Fig  4.8:  Effect  of  Reynolds  number  Rej}  on  separation  topology  (unsteady  cases). 

Shown  are  streamline  patterns  and  color  contours  of  streamwise  velocity  in  the  symmetry  plane 
(z=0,  on  top),  and  limiting  streamlines  and  color  contours  of  spanwise  vorticity  on  the  plate  (y-0, 
below).  Blue  indicates  reverse  flow. 
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4.3.4  Flow  Topology 


Foci  in  the  limiting  streamline  patterns  are  the  roots  of  so-called  “horn  vortices”,  as  described 
in  Sec.  1 .  The  topologies  observed  for  H  =  1 .0  and  Reynolds  numbers  greater  than  5,000  contain 
a  pair  of  foci  at  either  end  of  the  line  of  separation.  It  is  important  to  point  out  that  this  is  the  case 
for  both,  the  steady  and  unsteady  cases.  Furthermore,  for  all  cases  that  we  investigated  the 
streamline  pattern  in  the  symmetry  plane  displays  vortical  structures  with  reverse  flow  near  the 
surface. 

Insight  into  the  interaction  between  the  two  horn  vortices  associated  with  the  foci  on  the  flat 
plate  and  the  spanwise  vortex  associated  with  the  reverse  flow  in  the  symmetry  plane  can  be 
gained  from  3-D  streamlines  (Fig.  4.9a,b).  Fig.  4.9b  also  includes  a  visualization  of  the  vortex 
core,  which  we  extracted  using  the  method  of  Kenwright  and  Haimes  [1997J.  The  vortex  core 
extends  from  the  foci  upwards  and  towards  the  mid-span  plane.  Judging  from  the  visualizations 
in  Fig.  4.9.  the  horn  vortices  appear  to  connect  with  the  spanwise  vortex.  The  3-D  streamlines 
exhibit  a  strong  “swirling  motion”  induced  by  the  horn  vortices  and  also  illustrate  how  the  flow 
is  ingested  on  the  sides  and  then  ejected  near  the  top  of  the  separation  bubble  close  to  the 
symmetry  plane. 

When  analyzing  the  flowfield  in  cross-flow  planes  at  different  downstream  locations,  no 
streamwise  vortical  structures  could  be  identified  in  the  separated  region.  However,  as  shown  in 
Fig.  4.9c,  the  velocity  vectors  at  locations  downstream  of  the  separation  bubble  seem  to  indicate 
the  existence  of  a  pair  of  streamwise  vortices  close  to  the  surface. 
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(a)  Perspective  view  of  3-D  streamlines  and  color  contours  of  spanwise  vorticity  on  the  plate  obtained 
from  time-averaged  flow  data  for  Reo  ~  15,000,  W  =  6.4,  and  H  =  1.0  (blue  indicates  reverse  flow). 


(b)  Perspective  view  of  3-D  streamlines  and  color  contours  of  spanwise  vorticity  on  the  plate  for 
Reo  =  7,500 ,  W  =  12.8,  and  H  =  1.0  (steady;  blue  indicates  reverse  flow).  The  green  line  represents 
the  vortex  core  [Kenwright  &  Haimes  1997]. 


(c)  Cross-flow  visualization  for  Reo  =  7,500,  W  =  12.8,  and  H  =  1.0  downstream  of  separated  region 
at  x  =  0.6m  looking  downstream.  Velocity  vectors  and  color  contours  of  streamwise  vorticity  (red  - 
counter-clockwise  rotation,  blue  -  clockwise  rotation). 

Fig  4.9.  Flow  visualizations.  Cross-flow  planes  are  colored  with  isocontours  of 
streamwise  vorticity  (red  -  counter-clockwise  rotation,  blue  -  clockwise  rotation). 
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4.4  Results:  Unsteady  Separation  Phenomena 

Fig.  4.10  shows  isocontours  of  spanwise  vorticity  in  the  mid-span  plane,  for  the  time- 
averaged  flowfield  of  the  unsteady  cases  and  for  one  exemplary  steady  case  (Ren  =10,000).  For 
all  steady  cases,  the  separated  shear  layer  does  not  reattach  to  the  flat  plate.  For  the  unsteady 
cases,  however,  the  shear  layer  reattaches  in  the  mean.  Another  interesting  observation  is  that  as 
the  Reynolds  number  is  increased,  the  reattachment  point  moves  upstream,  resulting  in  a  gradual 
reduction  of  the  length  of  the  separation  bubble.  A  possible  explanation  is  provided  by  Fig.  4.1 1, 
which  shows  instantaneous  visualizations  of  the  spanwise  vorticity  in  the  mid-span  plane  for  the 
unsteady  cases.  The  separation  bubble  sheds  spanwise  vortical  structures,  which  break  up  into 
smaller  scale  structures,  leading  to  transition,  and  eventually  reattachment  in  the  mean.  These 
structures  are  likely  a  consequence  of  hydrodynamic  instabilities  of  the  flow  and  thus  highly 
dependent  on  the  flow  conditions.  As  the  Reynolds  number  is  increased  from  1 1,000  to  30,000, 
the  separated  boundary  layer  “rolls  up”  earlier  and  the  flow  transitions  farther  upstream,  resulting 
in  a  shortening  of  the  separation  bubble.  With  increasing  Reynolds  number,  the  size  of  the 
smallest  structures  is  reduced,  as  the  turbulent  energy  cascade  extends  over  a  larger  wave 
number  range.  In  addition,  the  wake  is  seen  to  spread  more  quickly  in  the  spanwise  direction. 
The  question  whether  the  primary  instability  is  convective  or  absolute  has  yet  to  be  addressed. 
Transition  not  only  results  in  an  earlier  closing  of  the  separation  bubble  but  also  affects  the 
instability  mechanisms  itself,  since  the  time-average  or  base  flow  is  altered.  In  summary,  the 
complicated  interplay  between  separation  and  transition  determines  the  separation  bubble 
characteristics. 

Wake  visualizations  are  provided  in  Fig.  4.12,  where  isosurfaces  of  the  Q  vortex 
identification  criterion  [Hunt  et  al.  1988]  are  shown.  Strong  coherent  structures  which  propagate 
in  downstream  direction  can  clearly  be  identified.  The  structures  have  a  strong  periodicity  in  the 
streamwise  direction,  even  for  Ren  =  30,000.  Similar  structures  were  observed  for  two- 
dimensional  separation  bubbles  (e.g.  Postl  2005).  These  large-scale  spanwise  coherent  structures 
are  mainly  responsible  for  the  transport  of  high-momentum  fluid  from  the  freestream  towards  the 
wall,  which  ultimately  results  in  reattachment  of  the  flow. 
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Fig  4.10:  Comparison  of  spanwise  vorticity  in  the  mid-span  plane  for  steady  and  unsteady 
separation  (for  the  latter  obtained  from  the  time-averaged  flow  data;  red  -  counter-clockwise 
rotation,  blue  -  clockwise  rotation). 
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Fig  4.11:  Comparison  of  instantaneous  spanwise  vorticity  in  the  mid-span  plane  for  cases  with 
unsteady  separation  (red  -  counter-clockwise  rotation ,  blue  -  clockwise  rotation). 
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Fig  4.12:  Visualization  of  vortical  structures  in  the  wake  of  the  separation  bubble.  Shown  are 
isosurfaces  of  the  Q  vortex  identification  criterion. 
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4.3.6  Comparison  to  Two-Dimensional  Separation  Bubbles 


Although  the  visualizations  of  spanwise  vorticity  in  the  symmetry  plane  provided  in  Fig.  4.10 
are  reminiscent  of  similar  visualizations  for  two-dimensional  separation  bubbles,  the  streamline 
patterns  in  the  symmetry  plane  (Figs.  4.4-4. 8)  are  evidence  that  no  closed  recirculation  region 
exists  for  the  steady  scenarios  that  we  investigated.  The  fluid  entrained  into  the  separation  bubble 
from  the  sides  is  ejected  at  the  top  of  the  bubble.  This  demonstrates  a  fundamental  difference 
between  two-  and  three-dimensional  separation  bubbles.  For  two-dimensional  bubbles,  the 
separated  boundary  layer  must  eventually  reattach,  thus  forming  a  closed  recirculation  region. 
Flowever,  in  three  dimensions,  the  attached  flow  downstream  of  the  separated  region  does  not 
require  a  closed  recirculation  in  the  symmetry  plane,  as  fluid  may  enter  the  reverse  flow  region 
from  the  sides.  For  the  steady  separation  scenarios  the  dividing  streamsurface  originating  from 
the  line  of  primary  separation  never  reattaches  (within  the  computational  domain).  Rather,  the 
attached  flow  downstream  of  the  separated  region  consists  of  fluid  that  is  being  entrained  from 
the  sides.  Flowever,  in  the  unsteady  cases  depicted  in  Fig.  4.10c-d,  the  shear  layer  of  the  time- 
averaged  flow  field  reattaches  just  as  for  a  two-dimensional  separation  bubble.  This  may  be  a 
consequence  of  the  large-scale,  unsteady  spanwise  coherent  structures,  as  discussed  in  Sec.  4.4. 

4.3.7  Comparison  to  Flow  over  Hemisphere-Cylinder 

The  wall  skin-friction  line  topology  suggested  by  Tobak  and  Peake  [1979]  for  the  separation 
bubble  at  the  nose,  or  “nose  bubble”  of  a  hemisphere-cylinder  at  low  to  intermediate  angles  of 
attack  is  depicted  in  Fig.  4.13.  It  features  a  saddle  point  of  separation  from  which  a  separation 
line  leads  to  a  pair  of  focal  points.  Two  saddle  points  and  a  node  of  attachment  in  the  symmetry 
plane  are  seen  in  the  rear  of  the  bubble.  The  foci  are  the  roots  of  two  “horn  vortices.  Lines  of 
secondary  separation  emerge  from  the  saddle  points  at  either  side  of  the  symmetry  plane  and 
extend  downstream. 

Our  results  for  the  unsteady  cases  (Fig.  4.8)  are  similar  to  the  limiting  streamline  topology 
proposed  by  Tobak  and  Peake  for  the  hemisphere-cylinder,  including  the  secondary  separation. 
For  the  steady  cases,  only  the  topologies  for  W  =  12.8  and  W  =  9.6  (Fig.  4.4c.d)  feature  the  same 
singular  points.  Both  cases,  however,  lack  any  secondary  separation. 

The  cross-flow  pattern  reported  by  Hsieh  and  Wang  [1996]  for  flow  over  a  hemisphere- 
cylinder  (Fig.  4.14)  features  two  streamwise  recirculation  regions  downstream  of  the  nose  bubble 
that  can  be  associated  with  the  cross-flow  component  of  the  approach  flow.  These  recirculation 
regions  are  similar  to  those  seen  in  our  simulations  of  3-D  separation  bubbles  (Fig.  4.9). 
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Fig.  4.13:  Skin-friction  line  topology  for  flow  over  hemisphere  cylinder  [Tohak  and  Peake  1979]. 
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5.  SEPARATION  BUBBLE:  WATER  TUNNEL  EXPERIMENTS 


5.1  Experimental  Setup 

In  parallel  with  the  numerical  investigations,  we  conducted  water  tunnel  experiments  for  the 
same  geometry.  We  considered  Reynolds  numbers  between  Reo=  5,000  and  Reo=l 5,000  based 
on  the  diameter  of  the  displacement  body  D.  A  schematic  of  the  experimental  setup  is  provided 
in  Fig.  5.1.  Boundary  layer  suction  was  applied  on  the  surface  of  the  displacement  body  to 
prevent  flow  separation  on  the  body  itself. 

SIDE  VIEW 


s 


FLOW 


> 


TOP  VIEW 


Duptocwitnt  lody  Ptot  plate 


Fig.  5.1:  Sketch  of  experimental  setup  showing  displacement  body ;  separation  bubble  and  stereo  P1V 
equipment  (not  to  scale). 


5.1.1  Exploratory  Experiments 

First,  a  preliminary  experimental  setup  was  developed  for  some  exploratory  experiments.  In 
order  to  speed  up  design  and  construction,  it  was  decided  to  construct  the  three-dimensional 
displacement  body  from  materials  that  were  already  available  in  the  lab  (Fig.  5.2).  This  allowed 
us  to  quickly  set  up  and  run  our  experiments  and  obtain  flow  visualizations  using  dye  injection  to 
get  a  first  glance  at  three-dimensional  separation  bubbles.  Based  on  these  observations,  we 
proceeded  with  the  design  of  an  advanced  and  accurate  displacement  body  with  improved 
boundary  layer  suction. 
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Fig.  5.2:  Preliminary  displacement  body.  Without  (left)  and  with  boundary  layer  suction  (right). 

5.1.2  Axisvmmetric  Displacement  Body 

The  displacement  body  was  built  from  a  fiberglass  shell.  A  negative  mold  was  generated 
from  a  positive  model,  which  had  been  fabricated  from  aluminum  (Fig.  5.3).  The  final  fiberglass 
displacement  body  was  obtained  using  the  negative  molds  (Fig.  5.4).  Both,  displacement  body 
and  flat  plate  were  painted  in  matte  black  to  reduce  glare  during  the  PIV  measurements.  To 
prevent  separation  from  the  displacement  body,  boundary  layer  suction  was  applied.  The  suction 
holes  had  a  diameter  of  1mm.  Gravity  forced  suction  was  generated  by  a  long  pipe  that  lead  to 
the  basement  of  the  building.  This  method  was  preferred  over  a  suction  pump  which  could 
introduce  undesired  disturbances.  The  removed  water  was  pumped  back  into  the  tunnel  to  ensure 
a  constant  static  pressure  in  the  test  section.  The  highest  suction  volume  flow  rate  was  10.5 
liter/min.  By  adjusting  the  length  of  the  suction  pipe  the  flow  rate  could  be  increased  or 
decreased.  Th^  adjustment  was  performed  whenever  the  flow  conditions  (e.g.  Reynolds  number) 
were  altered.  We  took  great  care  to  ensure  that  the  boundary  layer  suction  was  as  axisymmetric 
as  possible. 


Fig.  5.3:  Lathed  aluminum  positive. 
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Fig.  5.4:  Manufacturing  of  displacement  body  from  aluminum  positive  body ,  negative  fiberglass  mold  and 
positive  half  shell  (left).  The  final  displacement  body  with  boundary  layer  suction  holes  is  shown  on  the 
right. 
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Fig.  5.5:  Fluorescent  dye  injection  on  surface  of  modular  flat  plate. 


5.1.3  Flat  Plate 


In  addition  to  the  displacement  body,  a  flat  plate  model  was  manufactured.  The  flat  plate 
consists  of  a  leading  edge,  an  exchangeable  middle  section  with  dye  injection  holes  and  a  trailing 
edge.  The  surface  roughness  was  minimized  to  maintain  laminar  flow  over  the  entire  length  of 
the  flat  plate.  Interchangeable  middle-sections  allow  for  an  adjustment  of  the  plate’s 
configuration  and  length.  The  dimensions  of  the  flat  plate  were  chosen  such  that  the 
experimental  setup  could  be  varied.  The  dye  injection  holes  allow  us  to  visualize  the  flow  close 
to  the  wall.  The  dye  flow  rate  can  be  regulated  from  outside  the  water  tunnel.  A  typical  dye 
visualization  without  displacement  body  (Fig.  5.5)  clearly  indicates  that  the  flow  over  the  flat 
plate  is  laminar. 

5.2  Results  from  Dye  Flow  Visualization 

Before  we  carried  out  detailed  velocity  measurements  with  the  La  Vision  stereoscopic 
Particle  Image  Velocimetry  (PIV)  system,  we  employed  fluorescent  dye  for  flow  visualization. 
Using  dye  flow  visualizations  we  were  able  to  identify  the  reverse  flow  region  and  the  limiting 
streamlines  of  separation.  The  dimensions  of  the  separation  bubble  were  documented  for 
different  Reynolds  numbers.  Results  are  presented  for  a  Reynolds  number  range  of  Reo=5.000  to 
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Reo=l  5,000,  a  wall-normal  distance  of  the  displacement  body  from  the  flat  plate  of  H=l,  and  a 
downstream  distance  of  the  displacement  body  from  the  leading  edge  of  the  flat  plate  of  S  = 
2.5D.  For  a  Reynolds  number  of  Reo  =  5,000  the  bubble  was  found  to  be  steady.  For  Reo  = 
7,500  the  separation  bubble  was  almost  steady  (Fig.  5.6),  whereas  for  higher  Reynolds  numbers 
the  bubble  was  found  to  shed  vortical  structures. 


Fig.  5. 6:  Perspective  view  (left;  flow  direction  from  left  to  right)  and  view  from  the  back  (right;  flow 
direction  towards  the  observer)  for  ReD  =  7,500.  The  separation  bubble  was  visualized  using  fluorescent 
dye. 


Fig.  5.7:  Flow  visualization  without  (left)  and  with  (right)  boundary  layer  suction  for  Re»=7,500,  H-l 
and  S=2.5D. 

In  order  to  generate  a  sufficiently  large  pressure  gradient  on  the  flat  plate,  flow  separation 
from  the  displacement  body  had  to  be  prevented  by  applying  boundary  layer  suction.  Dye  flow 
visualization  was  utilized  to  optimize  the  suction  flow  rate  (Fig.  5.7).  Without  suction  the  flow 
separates  near  the  apex  of  the  displacement  body,  whereas  with  boundary  layer  suction  the  flow 
remains  attached  and  the  dye  lines  follow  the  contours  of  the  displacement  body. 

For  the  baseline  case  (Reo-7,500,  H=1  and  S=2.5D)  two  so-called  "horn  vortices"  are 
present  on  the  flat  plate.  The  topology  of  the  bubble  can  be  outlined  by  streamlines  in  the 
symmetry  plane  and  by  surface  skin  friction  lines.  The  best  flow  visualization  results  were 
obtained  when  the  fluorescent  dye  was  injected  through  four  holes  on  the  surface  of  the  flat  plate 
(Fig.  5.8).  The  limiting  streamlines  are  characterized  by  a  saddle  point  (S)  which  is  the  upstream 
separation  location  in  the  symmetry  plane  as  well  as  two  foci  (F)  which  are  situated  at  the  base  of 
the  horn  vortices.  To  capture  the  development  of  the  foci  in  our  visualizations,  only  the  outer 
streamlines  of  the  bubble  were  visualized  as  shown  on  the  right  in  Fig.  5.9. 
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Fig.  5.8:  Transient  development  of  separation  bubble  (after  pressure  gradient  is  applied  by  turning  on 
boundary  layer  suction)  for  the  baseline  case  (Ren=  7,500,  H=l,  S=2.5D). 


Fig.  5. 9:  Location  of  foci  in  three  dimensional  separation  bubble.  Flow  topology  sketch  (left)  by  Perry 
and  Chong  [1986]  and  dye  flow  visualization  (right)  for  Re  =  7,500,  H  =  1,  and  S  =  2.5D. 

For  high  enough  Reynolds  numbers  the  separated  flow  becomes  unstable,  resulting  in  the 
shedding  of  vortical  flow  structures  from  the  bubble.  With  increasing  Reynolds  number  the 
unsteady  flow  structures  appear  further  upstream,  leading  to  an  earlier  reattachment  of  the 
separated  flow.  This  trend  is  shown  in  Fig.  5.10.  We  also  found  that  by  bringing  the 
displacement  body  closer  to  the  flat  plate  the  shedding  frequency  could  be  increased.  In  addition, 
for  certain  conditions  we  observed  a  so-called  bubble  “breathing”  that  results  in  an  intermittent 
change  of  the  physical  size  of  the  separation  bubble  (Fig.  5.1 1). 
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Fig.  5.10:  Side  view  of  separation  bubbles  for  Reynolds  number  range  Re n  =  5,000  to  Ren  -  15,000  using 
fluorescent  dye.  As  the  Reynolds  number  is  increased  beyond  7,500,  the  bubble  begins  to  shed  vortical 
structures.  The  dominant  shedding  frequency  of  each  case  is  shown  in  the  lower  right  corner. 


Fig.  5.1 1 '  Sideview  of  the  separation  bubble  for  Re />= 5,000,  H  =  1,  and  S  =  2.5  D.  Dye  flow  visualizations 
showing  the  two  states  (steady  and  unsteady)  of  the  “ bubble  breathing'  cycle. 


5.3  PIV  Measurements 

For  all  Refolds  numbers  considered  (Reo  =  5,000  to  Reo  =  15,000)  we  measured  the 
instantaneous  velocity  components  in  the  mid-span  plane  and  parallel  spanwise  planes.  Since  the 
field  of  view  of  one  camera  is  limited  to  about  50mm,  several  sets  of  images  were  taken  and  later 
on  patched  together.  First,  instantaneous  velocity  measurements  were  obtained  at  a  number  of 
locations  in  downstream  direction,  and  then  the  averaged  velocity  field  over  the  entire  extent  of 
the  separation  bubble  was  constructed  from  the  instantaneous  data.  Finally,  velocity  profiles  at 
different  downstream  locations  were  extracted  from  the  time-averaged  PIV  data  (Fig.  5.12). 
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Normalized  with  the  freestream  velocity,  the  velocity  profiles  confirmed  the  different  local  flow 
behavior  which  was  already  observed  in  the  dye  flow  visualizations.  A  Blasius  velocity  profile 
was  measured  upstream  of  the  separation  bubble.  In  the  bubble,  inflectional  profiles  with  various 
degrees  of  reverse  flow  near  the  wall  were  obtained. 

Using  PIV,  we  also  measured  the  flow  around  the  displacement  body  to  once  more  confirm 
that  the  suction  strength  and  distribution  were  appropriate  (Fig.  5.13).  We  then  measured  the 
time-averaged  flow  in  the  mid-span  plane.  In  Fig.  5.14  the  results  from  these  measurements  for 
various  Reynolds  numbers  are  presented.  Streamlines  overlaid  with  streamwise  velocity 
contours,  as  well  as  velocity  profiles  at  different  downstream  locations  are  shown.  As  already 
observed  in  the  dye  flow  visualizations,  the  length  of  the  separation  bubble  decreases  with 
increasing  Reynolds  number.  In  addition,  the  time-averaged  PIV  measurements  also  indicate  that 
the  separation  bubble  height  decreases  with  increasing  Reynolds  number. 

For  a  Reynolds  number  range  of  Reo  =  7,500  to  Reo  =  15,000  we  observed  an  intermittency 
phenomenon.  Due  to  the  recirculating  fluid  in  the  separation  bubble,  the  size  of  the  bubble,  in 
particular  its  height,  increases.  The  separated  shear  layer  moves  away  from  the  wall  and  the 
reverse  flow  velocity  increases  to  a  point  where  the  bubble  starts  shedding  (Fig.  5.11).  The 
vortical  structures  associated  with  the  shedding  lead  to  flow  reattachment,  which  weakens  the 
reverse  flow  and  stops  the  shedding.  The  frequency  of  this  “breathing”  cycle  was  found  to  be  in 
the  order  of  fbreath  =  0.03  Hz  to  0.2  Hz  for  Reo  =  7,500  to  1 5,000. 
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Fig.  5.12:  PIV  data:  Time  averaged  velocity  profiles  at  different  downstream  locations  relative  to  the  flat 
plate  leading  edge  (normalized  with  free  stream  velocity ;  ReD  =5,000,  H=l,  S=2.5D). 
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Fig.  5. 13:  P1V  data:  Streamlines  showing  the  separated  boundary  layer  without  (top)  and  the  attached 
boundary  layer  with  boundary  layer  suction  (bottom). 
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Fig.  5.14:  PIV  data  for  Ren  =  5  000  to  Ren  =  15,000,  H  =  1  and  S  =  2.5D.  Streamlines  and  velocity 
contours  of  streamwise  velocity  in  the  mid-span  plane  (z=0),  as  well  as  velocity  profiles  at  different 
downstream  locations. 
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6.  SQUARE-DUCT  FLOW  SIMULATIONS 


Before  investigating  the  Stanford  diffuser  flow,  which  is  physically  more  complex  and  which 
features  a  turbulent  separation  bubble,  we  carried  out  simulations  of  a  square-duct  flow  using 
DNS,  RANS,  and  hybrid  RANS/LES.  The  purpose  of  these  simulations  was  to  determine  the 
grid  resolution  required  for  DNS  and  to  assess  RANS  and  hybrid  models  for  a  well  documented 
3-D  channel  flow.  For  these  simulations  the  Reynolds  number  based  on  bulk  velocity  and 
hydraulic  diameter  was  Reb=l  0,000. 


6.1  Computational  Grid 


The  hybrid  turbulence  model  simulations  were  carried  out  for  3  different  grid  resolutions  to 
determine  how  well  the  hybrid  models  adapt  to  changes  in  the  grid  resolution.  We  employed 
computational  grids  with  192*40x40  (coarse  grid,  307,200  cells),  1 92x54x54  (medium  grid, 
559,872  cells)  and  384x80x80  (fine  grid,  2,457,600  cells)  cells  (Fig.  6.1).  The  near-wall  grid 
resolution  in  wall  units  Ax+xAy+xAz+  was  82.1  x0. 941 X1 .88. ..43.6  (coarse  grid), 
41 .1  x0.941  xl  .88... 20. 6  (medium  grid),  and  20.5 x0. 941  xl  .88... 9. 87  (fine  grid).  In  addition,  as  a 
reference,  we  also  carried  out  a  direct  numerical  simulation  (DNS)  on  a  super  fine  grid  with 
384x160x160  cells  (9,830,400  cells.  Fig.  6.2)  where  we  resolved  all  scales  of  motion.  The  near 
wall  grid  line  spacing  in  wall  units  was  Ax+=10.3,  Ay+=0.314,  Az+=0.627...4.62.  The  grids  were 
generated  such  that  the  laminar  sublayer  was  resolved. 
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Fig.  6.1:  Cross-sectional  views  of  lower  left  corner  of  computational  domains  for  hybrid  RANS/LES 
simulations  of  square-duct  flow.  From  left  to  right:  Coarse ,  medium ,  and  fine  grid. 


Fig.  6.2.  Cross-sectional  view  of  computational  grid  for  DNS  of  square-duct  flow. 
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6.2  Domain  Length 

We  first  investigated  the  effect  of  the  domain  length.  Huser  and  Biringen  determined  a  zero 
streamwise  correlation  for  a  streamwise  separation  of  3.2  [Huser  and  Biringen  1993]. 
Accordingly,  a  domain  length  of  Lx  =  6.4  was  chosen.  Interestingly,  the  domain  length  in  the 
DNS  by  Gavrilakis  [1992]  was  Lx  =  5n  while  Madabhushi  and  Vanka  [Madabhushi  and  Vanka 
1991]  decided  on  a  domain  length  of  Lx  =  2n.  To  save  computer  time  we  performed  a  domain 
length  analysis  for  the  coarse  resolution  grid  using  FBR  [Johansen  et  al.  2004]  and  the  1998  k-co 
model  w  ith  EASM.  During  the  run  time  of  the  simulation  we  recorded  velocity  data  at  the  mid¬ 
channel  location.  Auto-correlations  of  the  recorded  u-velocity  component  are  showTi  in  Fig.  6.3. 
For  Lx  =  n  the  auto-correlation  has  multiple  peaks  which  can  be  associated  with  the  domain 
length  and  its  multiples.  As  the  domain  length  is  increased  by  factors  of  two  (with  the 
streamwise  grid  resolution  being  kept  constant  at  Ax  =  2tc/1  92)  “in-between”  peaks  disappear 
and  the  maxima  get  smaller  indicating  reduced  coherence  of  domain  length  related  waves.  For  Lx 
=  4n  a  weak  peak  can  still  be  discerned  for  At  *  10,  while  at  Lx  =  87c  the  auto-correlation  is 
almost  flat.  Isocontours  of  the  u-velocity  at  mid-span  (Fig.  6.4)  provide  visual  evidence  of  the 
large  structures. 

We  decided  that  a  domain  length  of  Lx  =  871  was  sufficient  for  our  coarse  grid  simulations. 
Out  of  practical  considerations  we  reduced  the  domain  length  to  4xc  for  the  medium  and  fine  grid 
and  to  2n  and  7t  for  the  super  fine  grid  to  keep  the  number  of  cells  and  the  computational  expense 
down  to  a  reasonable  limit. 


Fig.  6.3  Square-duct  flow  at  Re=10,000.  FBR  and  1998  k-co  model  with  EASM.  Auto-correlation  of  li¬ 
ve  loci  ty  disturbance  at  y  =  z  =  0.5  for  different  streamwise  domain  extents. 


Fig.  6.4.  Square-duct  flow  at  Re=10,000.  FBR  and  1998  k-co  model  with  EASM.  Isocontours  of  u-velocity 

(0  .  . .  1.5,  A  -  0.1)  at  mid-span. 
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6.3  Direct  Numerical  Simulation  (DNS) 


An  instantaneous  visualization  obtained  from  the  reference  DNS  for  Re=l 0,000  [Gross  and 
Fasel  2009a]  is  provided  in  Fig.  6.5.  Here,  the  Q  vortex  identification  criterion  by  Hunt  et  al. 
[1988]  was  employed  for  visualizing  vortical  structures.  The  near  wall  region  of  turbulence 
production  is  populated  by  worm-like  structures.  In  wall  bounded  flows  turbulent  streaks  have  an 
average  length  and  spanwise  spacing  in  wall  units  of  roughly  1000  and  100  [Cantwell  1981]. 


Fig.  6.5:  DNS  of  square-duct  flow.  Isosurfaces  of  Q-5. 


solid  lines 

dashed  lines  _ 

dotted  lines  _ 


Fig.  6.6:  DNS  of  square-duct  flow.  Isocontours  of  streamwise  velocity  in  cross-flow  plane  (left)  and  at 

mid  span,  z=0.5  (right). 


Fig.  6.7:  Auto-correlations  of  velocity  fluctuations  at  v=z=0.0493  (dotted  lines),  0.297  (dashed  lines), 

0.496  (solid lines). 


51 


Instantaneous  visualizations  of  the  streamwise  velocity  component  are  provided  in  Fig.  6.6. 
Also  indicated  are  3  points  which  were  situated  halfway  into  the  channel  at  y=z=0.0493,  0.297, 
0.496  where  time-dependent  velocity  data  was  extracted.  The  turbulent  structures  are  seen  to 
reach  all  the  way  to  the  center  of  the  channel.  Autocorrelations  of  the  velocity  data  obtained  at 
the  3  locations  indicated  in  Fig.  6.6  over  time  intervals  of  300  are  shown  in  Fig.  6.7.  For  the  short 
domain,  Ly  =  n,  peaks  are  seen  near  At  =  2.5  for  y  =  z  =  0.496  and  3.2  for  y  =  z  =  0.297.  With  a 
mean  flow  velocity  of  about  1 .32  near  the  center  of  the  channel  and  1 .25  near  y  =  z  =  0.297  the 
corresponding  downstream  wavelengths  are  3.3  and  4.  Therefore,  the  peak  at  At  =  2.5  is  likely 
related  to  the  domain  length.  Additional  confirmation  is  provided  by  the  fact  that  it  disappears 
for  the  longer  domain,  Lx  =  2 n.  A  peak  is  also  seen  near  At  =  4.8  for  both,  the  short  and  the  long 
domain.  This  peak  is  likely  related  to  the  domain  length  of  Ly  =  2n.  Although  our  domain  length 
study  showed  that  a  domain  length  of  271  was  insufficient,  out  of  practical  considerations  and  also 
because  Huser  et  al.  [1993]  used  a  domain  length  of  271  in  their  simulations,  we  did  not  further 
increase  the  downstream  extent  of  the  computational  domain  for  our  DNS. 


Fig.  6.8:  Energy  spectra  computed  from  velocity' fluctuations  for  short  (Lx—7t)  and  long  domain  (Lx=2n). 


We  also  computed  the  energy  spectra  (Fig.  6.8).  The  spectra  display  the  -5/3  slope  decay 
which  is  characteristic  of  the  inertial  subrange.  Distinct  peaks  are  visible  at  the  low  frequency 
end  of  the  spectrum.  Some  of  the  peaks  can  be  directly  related  to  the  peaks  in  the  autocorrelation, 
e.g.  for  the  short  domain  peaks  are  located  at  k  =  27tf  =  2n/2.5  =  2.5  and  2n/3.2  =  2.0.  Other 
peaks  may  likely  be  associated  with  energetic  coherent  flow  structures.  Since  the  mean  flow  is 
almost  grid  converged,  since  computations  with  longer  computational  domain  are  expensive,  and 
since  the  purpose  of  the  present  DNS  was  to  provide  mean  flow  reference  data  for  comparison 
with  RANS  and  hybrid  RANS/LES  results,  it  was  decided  that  the  long  domain  data  were 
sufficiently  accurate. 
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Fig.  6.9:  Temporal  and  streamwise  average  of  DNS  data  for  Lx=2n. 


Fig.  6.9  shows  temporal  and  streamwise  averages  of  the  flow.  As  a  result  of  the  finite 
duration  of  the  time-averaging  process  the  averages  are  not  fully  symmetric.  By  also  averaging 
over  all  octants  a  symmetric  picture  is  obtained.  When  considering  isocontours  of  the 
streamfunction,  T,  (v=5T/5z  and  w=-5T/3y)  counter-rotating  streamwise  vortices  emerge  that 
transport  high  momentum  fluid  from  the  center  into  the  comers.  This  secondary  flow  of  Prandtl’s 
second  kind  is  induced  by  the  Reynolds  stresses  and  typical  for  duct  flows  with  corners.  Almost 
identical  results  were  obtained  for  the  short  and  long  domain.  The  computed  bulk  velocities  of 
1 .03  for  the  short  and  long  domain  are  in  close  agreement  with  the  expected  result  of  1 .  The 
slight  mismatch  may  be  attributed  to  numerical  factors  such  as  accuracy  and/or  resolution,  the 
domain  length  or  the  empirical  relation  (we  used  Petukhov’s)  that  was  employed  for  obtaining 
the  magnitude  of  the  volume  force  that  drives  the  flow.  Velocity  profiles  in  wall  units  are 
compared  in  Fig.  6.10.  Our  DNS  data  are  in  good  agreement  with  other  simulation  data 
[Gavrilakis  1992,  Madabhushi  and  Vanka  1991,  Huser  and  Biringen  1993]. 
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Fig.  6.10:  Velocity  profiles  in  wall  units  at  wall  bisector. 


6.4  Reynolds-Averaged  Navier-Stokes  (RANS)  Calculations 

We  also  carried  out  steady  RANS  calculations  for  Re=  10,000  and  for  the  coarse,  medium, 
and  fine  grid  to  determine  a  suitable  underlying  RANS  model  for  the  hybrid  simulations. 
Because  hybrid  RANS/LES  reverts  to  RANS  in  the  coarse  grid  limit  the  accuracy  of  the  RANS 
model  is  important.  Isocontours  of  eddy  viscosity,  (It,  turbulence  kinetic  energy,  k,  and 
streamfunction,  'P,  taken  in  the  cross-flow  plane  are  shown  in  Fig.  6.11.  The  streamfunction 
isocontours  indicate  counter-rotating  streamwise  vortices  that  transport  high  momentum  fluid 
from  the  center  into  the  corners.  The  Explicit  Algebraic  Stress  Model  (EASM)  [Rumsey  and 
Gatski  2001]  is  needed  to  capture  this  secondary  flow.  With  the  low-Reynolds  number  terms  the 
near  wall  turbulence  kinetic  energy  is  increased.  A  comparison  with  the  DNS  data  is  provided  in 
Fig.  6.12.  Although  the  secondary  flow  is  captured  with  the  EASM  its  intensity  and  spatial 
characteristics  only  approximately  match  the  DNS  data. 


coarse 


medium 


fine 


1988  k-u  low- Re  k  EASM  1998  k-u)  k  EASM 
llr  k  'I'  tlr  k 


Fig.  6.11:  Square-duct  flow  at  Re= 10,000.  RANS  results .  Isocontours  of  eddy  viscosity  (pT  =  0  ..  .  40), 
turbulence  kinetic  energy  (k  =  0  .  .  .  0.02)  and  streamfunction  ( lower  left  quadrant ,  F  =  -0.002  .  .  . 
0.002)  in  cross-flow  plane  obtained  from  steady  RANS  calculations. 
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Fig.  6  12:  DNS.  Isocontours  of  instantaneous  u-velocity  (0<u<1.8,  Au=0.1).  Isocontours  of  u-velocity 
and  stream) unction  (-0.001<'F<0.001 ,  A'¥=0.0001)  for  R.  1  NS  results  obtained  with  1998  k-co  model  with 
Boussinesq  approximation  and  EASM. 

A  more  quantitative  comparison  becomes  possible  when  considering  profiles  taken  at  the 
wall  and  comer  bisector  and  comparing  with  our  DNS  results  and  the  DNS  data  by  Huser  and 
Biringen  [1993].  The  RANS  velocity  profiles  in  wall  units  (Fig.  6.13a)  all  overshoot  the  log- 
layer  law,  u  =  5  +  1/0.41  In  y+,  which  may  be  attributed  to  3-D  effects  resulting  from  the 
vicinity  of  the  neighboring  walls,  and  display  a  larger  log  layer  slope  than  the  DNS  data.  The 
mismatch  requires  further  attention  in  the  future.  Profiles  of  the  v-velocity  taken  at  the  comer 
bisector  allow  for  a  comparison  of  the  secondary  flow  amplitudes  (Fig.  6.13b).  Our  DNS  results 
agree  with  the  Huser  and  Biringen  [1993]  DNS  data  with  respect  to  the  maximum  amplitude, 
although  the  shapes  of  the  curves  are  slightly  different.  The  secondary  flow  is  under-predicted  in 
the  RANS  calculations  and  the  amplitude  distribution  does  not  match  the  DNS  data.  A 

comparison  of  the  Reynolds  normal  stresses,  u'u',  v'v',  and  w'w'  is  provided  in  Fig.  6.13c.  For 
our  DNS,  statistical  quantities  were  time-averaged  over  a  time  interval  of  400.  Agreement  of  our 
DNS  data  with  the  Huser  and  Biringen  [1993]  DNS  data  is  only  seen  for  the  streamwise  normal 

Reynolds  stress.  Discrepancies  in  the  v'v'  and  w'w'  distributions  necessitate  further 
examination.  The  RANS  data  match  the  DNS  data  only  qualitatively.  With  the  low-Reynolds 
number  terms,  higher  near  wall  normal  Reynolds  stresses  are  predicted.  Finally,  profiles  of  the 
Reynolds  shear  stress  and  the  turbulence  kinetic  energy  are  shown  in  Fig.  6.13d.  Surprisingly, 
the  Reynolds  shear  stress  profiles  from  the  DNS  and  RANS  are  quite  close.  A  very  good  match 
between  the  two  DNS  data  sets  is  seen  for  the  turbulence  kinetic  energy.  Especially  near  the  wall 
the  turbulence  kinetic  energy  distribution  obtained  from  the  RANS  calculation  with  low- 
Reynolds  number  terms  agrees  reasonably  well  with  the  DNS  data.  Without  low-Reynolds 
number  terms,  the  near  wall  turbulence  kinetic  energy  peak  is,  as  expected,  under-predicted. 
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Fig.  6.13:  Square-duct  flow  at  Re=10,000.  a)  Wall  bisector  profiles  of  streamwise  velocity  in  wall  units 
b)  Corner  bisector  profiles  of  v-velocity  component,  c)  and  d)  Wall  bisector  profiles  of  Reynolds  stresses 
and  turbulence  kinetic  energy.  Black  lines:  u  from  Fig.  6,  run  B.  — v  from  Fig.  7b,  run  B,  u'u'  and  u'v' 
from  Fig.  llc/a,  vV  and  w'w'  from  Fig.  13a/b  in  Huser  and  Biringen  [1993];  Green  lines:  our  DNS; 
l  'van  lines :  RANS,  1988  k-co  low-Re  and  EASM;  Blue  lines  RANS,  1998  k-co  and  EASM. 

6.5  Hybrid  RANS/LES  Simulations 

We  then  computed  the  same  flow  using  two  hybrid  turbulence  models  that  we  determined 
successful  in  an  earlier  numerical  campaign  [Gross  and  Fasel  2008b,  2009a, b].  We  employed 
both,  a  modified  version  of  our  flow  simulation  methodology  (FSM)  [Speziale  1997,  Fasel  et  al. 
2002]  where  the  model  contribution  is  computed  from  the  von  Karman  energy  spectrum  and  the 
filter-based  RANS  approach  [Johansen  et  al.  2004]  with  filter-width  set  equal  to  the  local  grid 
line  spacing.  An  approach  similar  to  the  one  by  Batten  et  al.  [2004]  was  employed  for  “injecting'5 
or  “seeding”  statistically  random  turbulent  motion  in  regions  with  vanishing  model  contribution 
[Gross  and  Fasel  2009b].  As  underlying  turbulence  models  we  employed  the  1988  k-co  model 
with  low-Reynolds  number  terms  and  the  1998  k-co  model,  both  with  EASM.  Time-averages 
were  computed  over  time-intervals  of  400,  except  for  the  coarse  grid  FSM  with  1988  k-co  model 
and  low-Reynolds  number  terms  and  the  medium  grid  FSM  with  1998  k-co  model  for  which 
time-averages  were  computed  over  time-intervals  of  360  and  300,  respectively. 

Velocity  fluctuations  were  computed  with  respect  to  the  unprocessed  “raw”  time-averages, 
meaning  that  the  time-averages  were  not  averaged  in  the  streamwise  direction  and  over  the  8 
octants.  The  reasoning  behind  this  approach  was  that  any  imbalances  or  inaccuracies  in  the  code 
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that  may  result  in  an  asymmetry  in  the  solution,  which  would  be  captured  in  the  time-averages, 
would  result  in  a  systematic  error  if  the  time-averages  were  averaged  in  the  streamwise  direction 
or  over  the  8  octants.  For  the  cases  where  the  computed  bulk  velocities  were  close  to  the  DNS 
value  of  1 .03  statistical  quantities  were  computed  over  a  time  interval  of  400. 
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Fig.  6.14:  Square-duct  flow  at  Re=10,000.  Hybrid  RANS/LES  and  DNS.  Isosurfaces  of  Q  =  0.1,  and 
isocontours  of  streamj unction,  *¥  =  -0.002  .  .  .  0.002 ,  resolved  turbulence  kinetic  energy,  kr-  0  .  .  .  0.02, 
unresolved  turbulence  kinetic  energy,  ku  =  0  .  .  .  0.02,  model  contribution,  f  =  0  ...  1,  and  unresolved 
eddy  viscosity,  juju  —  0  ...  20. 


An  overview  of  the  different  cases  is  provided  in  Fig.  6.14.  When  considering  the  isosurfaces 
of  the  vortex  identification  criterion,  Q  =  0.1,  which  allows  for  an  identification  of  vortical  flow 
structures,  the  most  small-scale  structures  are,  as  expected,  observed  for  the  DNS  on  the  super 
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fine  grid.  As  the  grid  resolution  is  reduced  more  and  more  of  the  resolved  turbulence  kinetic 
energy,  kr,  is  absorbed  in  the  unresolved  turbulence  kinetic  energy,  ku,  and  less  and  less  small- 
scale  flow  structures  are  resolved.  Isocontours  of  the  streamfunction,  'F,  taken  in  a  cross-flow 
plane  illustrate  the  secondary  flow.  Differences  in  the  predicted  secondary  flow  distributions  are 
indicative  of  differences  in  the  particular  Reynolds  stress  distributions  which  drive  the  secondary 
flow  (Fig.  6.14).  Huser  et  al.  [1994]  compared  results  obtained  from  a  RANS  calculation  based 
on  Speziale’s  non-linear  k-8  model  with  DNS  data.  They  found  significant  differences  between 
the  model  prediction  and  the  DNS  results,  such  as  an  underprediction  of  the  intensity  of  the 
secondary  flow.  The  model  contribution,  f,  provides  a  measure  for  how  much  turbulence 
modeling  is  afforded  by  the  hybrid  RANS/LES  models.  With  increasing  grid  resolution  the 
model  contribution  is  reduced.  With  FBR  and  1988  k-co  model  with  low-Reynolds  number  terms 
the  model  contribution  is  1  at  the  highest  grid  resolution  and  the  unresolved  eddy  viscosity  is 
close  to  0.  This  inconsistency  requires  further  investigation.  In  general,  a  lower  model 
contribution,  f,  coincides  with  a  smaller  unresolved  eddy  viscosity,  piu,  and  more  resolved  flow 
structures.  In  all  instances,  the  model  contribution  is  1  near  the  wall  which  is  desirable  as  modem 
RANS  models  capture  wall  bounded  turbulent  flows  with  high  fidelity.  The  model  contribution 
is  seen  to  be  slightly  larger  for  the  1998  k-co  model  than  for  the  1988  k-co  model  with  low- 
Reynolds  number  terms.  Also,  the  model  contribution  is  lower  for  the  FSM  than  for  the  FBR 
which  again  manifests  itself  in  a  larger  number  of  resolved  flow  structures.  However,  for  the 
FBR  model  the  constant  C3  could  be  adjusted  for  reducing  the  model  contribution.  The 
unresolved  eddy  viscosity,  p.xu,  even  on  the  coarse  grid,  is  always  less  than  20  which  indicates 
that  the  turbulence  model  contribution  in  the  present  simulations  is  relatively  lowr. 
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Tab.  6.1 :  Square-duct  flow  at  Re= 10,000.  Computed  bulk  velocities.  *-  Discontinued,  approximate 

results;  t-  Steady  solutions. 

A  summary  of  the  computed  bulk  velocities  is  provided  in  Tab.  6.1.  The  DNS  result  is 
highlighted  in  green  as  it  provides  the  reference  for  the  hybrid  simulations.  Since  we  use  the 
same  code  for  the  hybrid  simulations,  the  DNS  result  should  be  approached  in  the  fine  grid  limit. 
The  RANS  results  are  highlighted  in  blue  as  they  provide  the  low-grid  resolution  reference  for 
the  hybrid  simulations.  As  the  grid  resolution  is  reduced  the  RANS  result  should  be  approached. 
The  bulk  velocities  obtained  from  our  hybrid  simulations  with  FBR,  1998  k-co  model,  and 
EASM  approach  the  RANS  result  in  the  coarse  grid  limit  and  the  DNS  result  in  the  fine  grid 
limit.  With  FBR,  1988  k-co  model.  low-Reynolds  number  terms,  and  EASM  a  very  good  bulk 
velocity  prediction  is  already  obtained  for  the  medium  resolution  grid,  however,  for  the  fine  grid, 
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as  was  seen  earlier  (Fig.  6.14)  the  unresolved  eddy  viscosity  goes  to  zero  and  the  predicted  bulk 
velocity  is  too  high.  With  1988/98  k-co  model,  Boussinesq  approximation,  and  low-Reynolds 
number  terms  steady  solutions  were  obtained  for  the  coarse  grid.  Medium  grid  resolution 
simulations  were  not  attempted.  The  fine  grid  simulation  with  1998  k-co  model  was  unsteady. 
With  the  present  version  of  the  FSM  the  predicted  bulk  velocities  are  slightly  too  high  and  the 
model  contribution  is  too  low. 


Fig.  6.15:  Square-duct  flow  at  Re=l  0,000.  Isocontours  of  streamwise  vorticity,  cox  =  ~10.  ..10,  at z  =  0 
obtained  from  FBR  with  1998  k-co  model  and  EASM  and  DNS. 

Isocontours  of  the  streamwise  wall  vorticity, 


co. 


dw  dv 
dy  dz  ’ 


(6.1) 


provide  evidence  of  near  wall  streamwise  vortical  structures  (Fig.  6.15).  Again  it  can  be  seen 
how  the  number  and  intensity  of  the  flow  structures  is  reduced  with  decreasing  grid  resolution. 
Velocity  profiles  and  profiles  of  turbulence  kinetic  energy  (unresolved,  resolved,  and  total)  for 
the  FBR  are  shown  in  Fig.  6.16.  The  medium  and  fine  grid  velocity  profiles  obtained  with  both. 
1988  k-co  model  with  low-Reynolds  number  terms  and  1998  k-co  model,  almost  match  the  DNS 
profile  in  the  log-layer.  The  coarse  grid  velocity  profiles  approach  or  overshoot  the  RANS  result 
near  the  center  of  the  channel.  The  total  turbulence  kinetic  energy,  k,  was  computed  by 
summation  of  the  time-averaged  unresolved  turbulence  kinetic  energy,  ku,  which  attains  its 
maximum  near  the  wall,  and  the  time-averaged  resolved  turbulence  kinetic  energy,  kr,  which  was 
computed  from  the  unsteady  velocity  fluctuations  and  which  attains  its  maximum  further  away 
from  the  wall  and  more  towards  the  center  of  the  channel.  The  total  turbulence  kinetic  energy  is 
higher  for  the  1988  k-co  model  with  low-Reynolds  terms  than  for  the  1998  k-co  model.  With  1998 
k-co  model  the  total  turbulence  kinetic  energy  shows  the  same  near  wall  magnitude  as  for  the 
DNS.  As  the  grid  resolution  is  reduced  the  location  of  the  peak  of  the  total  turbulence  kinetic 
energy  moves  away  from  the  wall  as  the  unresolved  turbulence  kinetic  energy  distribution  begins 
to  approach  the  RANS  result. 
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a) 

FBR,  1988  k- u  low- Re  k  EASM 


b ) 


FBR,  1998  k-oj  k  EASM 


Fig.  6.16:  Square-duct  flow  at  Re= 10,000.  Wall  bisector  profiles  of  a)  streamwise  velocity  in  wall  units 
and  b)  total  (red),  resolved  ( >rang  )  and  unresolved  (magenta)  turbulence  kinetic  energy'  obtained  from 
FBR  simulations..  Green  lines  DNS;  ('van  lines  RANS,  1988  k-m  low-Re  and  EASM;  Blue  lines:  RANS, 
1998  k-co  and  EASM.  Red  lines:  hybrid  models  as  indicated.  Dotted  lines:  Coarse  grid;  Dashed  lines: 
Medium  grid;  Solid  lines.  Fine  grid. 
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n) 

FSM,  1988  k-oj  low-Re  k  EASM 


b) 


FSM.  1998  k-u  k  EASM 


Fig.  6.17;  Square-duct  flow'  at  Re= 10,000  Wall  bisector  profiles  of  a)  streamwise  velocity  in  wall  units 
and  b)  total  (red),  resolved  (orange)  and  unresolved  (magenta)  turbulence  kinetic  energy >  obtained  from 
FSM  simulations .  Green  lines  DNS;  Cyan  lines  RANS '  1988  k-co  low-Re  and  EASM;  Blue  lines:  RANS, 
1998  k-co  and  EASM.  Red  lines  hybrid  models  as  indicated.  Dotted  lines:  Coarse  grid,  Dashed  lines: 
Medium  grid;  Solid  lines :  Fine  grid. 

Similar  graphs  for  the  FSM  are  provided  in  Fig.  6.17.  Turbulence  kinetic  energy  distributions 
are  not  shown  for  the  FSM  with  1988  k-co  model  and  low-Reynolds  terms.  Because  the  predicted 
bulk  velocities  were  too  high  we  decided  against  computing  statistical  quantities.  The  velocity 
profiles  in  wall  units  obtained  from  the  FSM  are  closer  to  the  RANS  than  to  the  DNS  result  as 
the  grid  resolution  is  increased  from  coarse  to  medium.  Turbulence  kinetic  energy  distributions 
are  only  available  for  the  coarse  grid  FSM  with  1998  k-co  model.  Compared  with  the  FBR  result, 
the  resolved  turbulence  kinetic  energy  and  accordingly  also  the  total  turbulence  kinetic  energy 
are  increased. 
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a)  b) 

FBR,  1988  k-uj  low-Re  &  EASM 


Fig.  6  18:  Square-duct  flow  at  Re— 10,000.  a)  Auto-correlations  and  b)  frequency  spectra  of  u-velocity 
perturbation  at  y  =  z  =  0.5.  Green  lines:  DAS;  hybrid  models  as  indicated:  Red  lines  Coarse  grid; 
Magenta  lines  Medium  grid;  Orange  lines  Fine  grid.  Black  lines.  -5/3  slope  of  inertial  sub-range. 


Auto-correlations  and  frequency  spectra  of  the  u-velocity  perturbation,  u’=u-U,  at  the  mid¬ 
channel  location  are  shown  in  Fig.  6.18.  The  data  appears  sufficiently  de-correlated  in  time, 
leading  to  the  conclusion  that  the  streamwise  domain  extent  was  likely  sufficient.  However,  the 
original  intent  was  to  identify  extrema  in  the  auto-correlations  and  correlate  them  with  dominant 
wavelengths  in  the  flow  and  then  to  investigate  how  well  those  waves  are  represented  by  the 
different  hybrid  RANS/LES  approaches  for  different  grid  resolutions.  Spectra  are  shown  in  Fig. 
6.18b.  The  low  frequency  end  of  the  spectra  is  not  well  resolved  because  of  the  limitation  of  the 
time  intervals  over  which  time-dependent  data  was  recorded.  Nevertheless,  Fig.  6.18b  indicates 
that  the  inertial  sub-range  is  captured  by  the  DNS  and  partially  captured  by  the  medium  and  fine 
grid  hybrid  simulations.  Fig.  6.18b  also  shows  how  the  wave  number  cutoff  is  shifted  to  lower 
wave-numbers  as  the  grid  resolution  is  reduced. 
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6.6  Rectangular  Duct  Flow 

As  a  next  step  towards  simulations  of  the  entire  diffuser  we  then  computed  the  approach 
channel  flow  of  the  Stanford  diffuser  experiment  [Cherry  et  al.  2008].  The  approach  channel 
flow  Reynolds  number  based  on  channel  height  was  10,000  (or  Re=l 5,381  based  on  hydraulic 
diameter).  The  computational  domain  (Fig.  6.19)  had  dimensions  (lengthx  heightx  width)  of 
6.66x1x3.33  and  a  resolution  of  64x96x96  cells.  The  near  wall  grid  line  spacing  in  wall  units 
was  Ax+=62,  Ay+=0.89,  and  Az+=1.8...33  at  the  “long”  wall  (y=0)  and  Ax+=62,  Ay+=1.78...7.1, 
and  Az+=0.89  at  the  “short”  wall  (z=0).  Again,  we  first  obtained  RANS  solutions  with  the  1998 
k-co  model  (Fig.  6.20).  In  agreement  with  the  experiment  a  secondary  flow  is  captured  with  the 
EASM  (Fig.  6.21).  The  computed  bulk  velocities  are  1.05  (Boussinesq  approximation)  and  1.08 
(EASM). 


Fig.  6. 1 9:  Computational  grid for  RANS  and  hybrid  simulations  of  approach  channel  flow. 
time-av.  u-velocity  streamf unction 

1998  k-(0 


1998  k-co,  EASM 


Fig.  6.20:  RANS  results  obtained  with  1998  k-co  model  Isocontours  of  streamwise  velocity  (0<u<L8, 
Au=0.1)  and  stream  function  (-0.00 1  <¥<0.001 ,  A¥=0.0001).  Top:  Boussinesq  approximation,  bottom: 
EASM. 


Instantaneous  flow  visualizations  obtained  from  hybrid  simulations  with  FSM  and  filter- 
based  RANS  (FBR)  are  shown  in  Fig.  6.22.  We  employed  both  approaches  with  and  without 
“turbulence  injection”  in  areas  of  vanishing  model  contribution.  Slightly  fewer  structures  are 
seen  for  the  FBR  again,  indicating  a  larger  model  contribution.  No  difference  is  seen  between  the 
instantaneous  flow  fields  with  and  without  “turbulence  injection”. 
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Fig.  6.21:  Velocity  profiles  at  wall  bisectors.  Circles:  Experiment  by  Cherry  et  al.  [2008], 


FSM  filter-based  RANS 


Fig.  6.22:  Hybrid  simulations  of  approach  channel  flow  without  (top)  and  with  “ turbulence  injection  ” 
(bottom).  Isocontours  of  Q=0.1.  Left:  FSM  and  right:  FBR 


Fig.  6.23:  Approach  channel  flow  bulk  velocities. 

Histories  of  the  computed  bulk  velocities  are  shown  in  Fig.  6.23.  The  bulk  velocities  are 
generally  higher  for  the  FSM  when  compared  to  the  filter-based  RANS.  “Turbulence  injection” 
lowers  the  bulk  velocities  possibly  as  the  result  of  more  and  strengthened  small  scale  structures. 
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7.  DIFFUSER  FLOW  SIMULATIONS 


We  considered  our  hybrid  RANS/LES  duct  flow  simulations  successful  enough  to  begin 
simulations  of  a  geometrically  and  physically  more  complex  flow,  the  Stanford  asymmetric 
diffuser  for  which  measurements  indicate  a  turbulent  separated  flow  region  [Cherry  et  al.  2008]. 
The  diffuser  inflow  bulk  Reynolds  number  based  on  hydraulic  diameter  is  Re=l  5,380. 

7.1  Computational  Grid 

The  computational  domain  for  our  simulations  was  broken  up  into  two  sub-domains,  one  for 
the  approach  channel  flow  and  one  for  the  diffuser  (Fig.  7.1).  The  approach  channel  flow  sub- 
domain  had  a  length  xheightxwidth  ratio  of  6.66x1x3.33  with  92x76x96  cells.  The  near  wall  grid 
line  spacing  in  wall  units  was  Ax+  =  42.8,  Ay+=0.556,  and  2.02  <  Az+  <  36.3  at  y=0  and  Ax+  = 
42.8,  1.07  <  Ay+  <  13.6,  and  Az+  =  1.01  at  z  =  0.  The  diffuser  opens  up  from  a  heightx  width 
ratio  of  1  x3.33  at  x  =  0  to  4X4  at  x  =  15.  The  inflow  of  the  diffuser  sub-domain  was  located  at  x 
=  -0.579  to  not  fully  suppress  potential  upstream  effects  of  the  diffuser  flow  on  the  approach 
channel  flow.  The  outflow  was  located  at  x  =  45.  The  diffuser  grid  had  268x76x96  cells.  The 
near  wall  grid  line  spacing  in  wall  units  at  the  diffuser  exit  (Re=8,325)  was  Ax+  =  18.5,  Ay+  = 
0.473,  and  0.548  <  Az+  <  9.87  at  y  =  0  and  Ax+  =  18.5,  0.946  <  Ay+  <  12.0,  and  Az+  =  0.274  at  z 
=  0.  The  total  number  of  cells  was  2.7  million. 


4x4 


Fig.  7.1:  Computational  grid  for  diffuser  flow  simulations. 

7.2  Simulation  Results 

We  carried  out  several  hybrid  simulations  [Gross  and  Fasel  2008b,  2010].  According  to 
Margolin  and  Rider  [2002]  certain  second-order-accurate  upwind  schemes  have  similar 
properties  as  standard  LES  sub-grid  models.  This  motivated  us  to  perform  a  simulation  with 
second-order-accurate  TVD  scheme  [Yee  1987]  for  the  convective  terms  and  a  second-order- 
accurate  discretization  for  the  viscous  terms.  In  addition,  we  performed  one  simulation  with  the 
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same  discretization  as  for  the  hybrid  simulations  (ninth-order-accurate  discretization  for  the 
convective  terms  and  fourth-order-accurate  discretization  for  the  viscous  terms)  but  without 
turbulence  model.  The  time  step  for  all  4  simulations  was  At  =  0.01.  The  simulations  without 
turbulence  model  were  advanced  in  time  over  12000  time  steps,  which  corresponds  to 
approximately  7  diffuser  flow-through  times.  For  the  FBR  hybrid  simulation  with  1988  k-co 
model,  low-Reynolds  number  terms,  and  EASM  time-averaging  was  initiated  after  28000  time 
steps  and  stopped  after  36000  time  steps  (approximately  5  flow-through  times).  For  the  FBR 
hybrid  simulation  with  1998  k-w  model  and  EASM  the  time-averaging  was  initiated  after  12000 
time  steps  and  terminated  after  36000  time  steps  (approximately  14  flow-through  times). 


Fig.  7.2:  Diffuser  approach  flow  bulk  velocities. 


The  time-histories  of  the  approach  flow  bulk  velocity  are  depicted  in  Fig.  7.2.  With  the 
second-order-accurate  TVD  scheme  all  unsteady  flow  structures  in  the  approach  flow  are 
dampened  out  hinting  at  excessive  numerical  diffusion.  This  loss  of  unsteady  flow  structures 
reduces  the  amount  of  wall-normal  mixing  which  increases  the  bulk  velocity.  Without  correction 
of  the  channel  flow  volume  force  the  bulk  velocity  increases  beyond  1 .8.  The  resulting  diffuser 
inflow  Reynolds  number  is  larger  than  in  the  experiment.  With  the  higher-order-accurate  scheme 
and  without  turbulence  model  the  approach  flow  bulk  velocity  asymptotically  approaches  a  value 
in  excess  of  1.3.  With  the  hybrid  RANS/LES  models  the  approach  flow  bulk  velocity  stabilizes 
near  1.04  which,  considering  the  value  of  1.03  for  the  square-duct  flow  DNS.  is  an  acceptable 
value. 
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Fig.  7.3.  Isosurfaces  of  Q-0.01  colored  by  u-velocity  and  isocontours  of  u-velocity  (0...1.5)  at  z= 1.665 

Instantaneous  flow  visualizations  using  the  Q-criterion  [Hunt  1988]  are  shown  in  Figs.  7.3 
and  7.4.  In  all  instances  small-scale  unsteady  flow  structures  are  seen  in  the  diffuser  (0  <  x  <  15) 
even  for  the  diffusive  second-order-accurate  scheme  where  all  unsteadiness  in  the  approach 
channel  flow  is  suppressed.  This  hints  at  a  strong  hydrodynamic  instability  mechanism  that 
generates  unsteady  flow  structures.  Few  large  structures  are  seen  for  the  diffusive  second-order- 
accurate  scheme  while  many  small  structures  are  seen  for  the  ninth-order-accurate  scheme 
without  turbulence  modeling. 
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Fig.  7.4.  Diffuser  at  Re=l 5,380.  Isosurfaces  of  Q  =  /.  Side  views  and  top  down  views. 


Spanwise  views  of  the  time-averaged  streamwise  velocity  (Fig.  7.5)  provide  information 
regarding  the  shape  and  extent  of  the  separated  flow  region.  In  the  experiments,  volume  data 
were  obtained  using  magnetic  resonance  velocimetry  [Cherry  et  al.  2008].  In  the  spanwise  plane 
z=l  .665  the  flow  appears  to  separate  near  x  »  7.5  and  not  at  the  diffuser  inlet  comer  at  x  =  0.  A 
similar  behavior  is  observed  for  the  FBR  with  1988  k-co  model,  low-Reynolds  number  terms,  and 
EASM,  although  no  conclusive  statement  can  be  made  yet  because  the  quality  of  the  time- 
averages  is  poor.  Because  the  diffuser  geometry  is  asymmetric,  data  can  only  be  averaged  in  time 
and  time-averages  take  a  long  time  to  converge.  For  the  FBR  with  1998  k-co  model  the  flow 
appears  to  separate  at  the  comer  where  the  diffuser  begins  to  open  up,  (x,  z)  =  (0,  1 .665),  and  to 
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reattach  near  (x,  z)  =  (7.5,  1 .665).  Also  noteworthy  is  that  the  approach  channel  flow  bulk 
velocity  appears  to  be  smaller  in  the  experiment. 
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Fig.  7.5:  Diffuser  at  Re~ 15,380.  Isocontours  of  time-averaged  streamwise  velocity  at  z  =  1.665. 

Measurements  by  Cherry  at  al.  [200%]. 
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Fig.  7.6:  Diffuser  at  Re=l 5,380.  Isocontours  of  instantaneous  unresolved  turbulence  kinetic  energy , 
ku=0  .  .  .  0.02,  model  contribution^  =  0  .  .  1,  and  unresolved  eddy  viscosity ,  /uTu  =  0  ...  10,  at  z  =  1.665. 


Instantaneous  visualizations  in  the  spanwise  plane,  z  =  1 .665,  for  the  two  hybrid  RANS/LES 
simulations  are  shown  in  Fig.  7.6.  Unresolved  turbulence  kinetic  energy,  ku,  is  produced  in  the 
separated  shear  layer  and  near  the  walls.  Both,  model  contribution,  f,  and  unresolved  eddy 
viscosity,  piu-  are  larger  for  the  1998  k-co  model  than  for  the  1988  k-co  model  with  low-Reynolds 
number  terms. 
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Fig  7. 7:  Diffuser  at  Re=l 5,380.  Cross-sectional  views  of  streamwise  velocity >,  U-—0.2.  .  1. 

Cross-sectional  views  of  streamwise  velocity  as  shown  in  Fig.  7.7  provide  information  about 
the  shape  and  extent  of  the  separated  flow  region.  In  the  experiment  at  x=5,  the  flow  was 
separated  from  the  top  right-hand-side  comer.  At  x=10  and  x=15  the  flow  was  separated  from 
the  entire  top  wall  with  the  velocity  maximum  being  slightly  off  to  the  left  from  mid-span.  At 
x=20  the  flow  was  fully  attached.  Although  the  quality  of  the  time-averages  obtained  from  our 
hybrid  simulations  is  insufficient  to  allow  for  any  final  conclusions,  the  amount  of  flow 
separation  appears  under-predicted  compared  to  the  experiment. 
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8.  SUMMARY  AND  CONCLUSIONS 


8.1  3-D  Separation  Bubble  on  Flat  Plate:  Simulations  and  Experiments 

We  studied  steady  and  unsteady  3-D  separation  bubbles  on  a  flat  plate  using  both  direct 
numerical  simulations  (DNS)  and  water  tunnel  experiments.  Separation  was  induced  by  a  three- 
dimensional  pressure  distribution  that  was  generated  by  an  axisymmetric  displacement  body 
which  was  placed  above  the  flat  plate.  Fundamental  features  that  are  documented  in  the  literature 
for  3-D  separated  flows  in  general  as  well  as  for  the  flow  over  a  hemisphere-cylinder  in 
particular  could  be  reproduced  in  our  simulations  and  experiments.  We  investigated  the  effects  of 
changes  in  the  Reynolds  number  and  the  non-dimensional  pressure  distribution  on  the  topology 
of  the  separation.  The  closed  separation  topologies  that  we  identified  all  showed  the  same 
characteristic  pattern  of  the  skin  friction  lines  at  the  upstream  end  of  the  separation  bubble;  i.e., 
the  line  of  separation  connects  a  saddle  point  of  separation  with  a  focus  on  either  side  of  the 
symmetry  plane.  The  foci  are  the  roots  of  so-called  “horn-vortices”.  We  also  showed  that  the 
topology  at  the  downstream  end  of  the  separation  bubble  varies  considerably  for  different  flow 
conditions.  In  addition  to  the  closed  separation  patterns,  we  also  found  an  open  separation 
topology  that  had  not  been  documented  before.  The  separation  patterns  that  we  identified  for  the 
unsteady  cases  exhibit  a  topology,  which  is  equivalent  to  what  is  documented  in  the  literature  for 
the  separated  flow  over  a  hemisphere-cylinder  at  low  to  intermediate  angles  of  attack.  In 
particular,  we  could  detect  lines  of  secondary  separation.  The  fact  that  the  characteristics  of  the 
separated  flow  over  a  hemisphere-cylinder  could  not  always  be  reproduced  in  its  entirety  in  our 
flat  plate  simulations,  suggests  that  future  research  needs  to  address  the  effect  of  surface 
curvature.  Our  results  also  suggest  that  by  adjusting  the  wall  pressure  distribution  it  is  possible  to 
at  least  qualitatively  duplicate  the  flow  behavior  at  a  different  Reynolds  number.  Our  simulation 
results  are  in  qualitative  agreement  with  the  accompanying  water  tunnel  experiments.  In 
particular,  the  topology  of  the  bubble  and  the  onset  of  bubble  shedding  at  larger  Reynolds 
numbers  could  be  confirmed.  We  also  noticed  a  relationship  between  Reynolds  number,  vortex 
shedding  and  the  so-called  “bubble  breathing”,  as  observed  by  Dogval  and  Kozlov  [1994]  for 
two-dimensional  separation  bubbles.  The  experiments  thus  provided  guidance  and  validation  for 
our  simulations. 


8.2  Simulations  of  Square-Duct  Flow  and  Asymmetric  Diffuser 

We  computed  a  square-duct  flow  at  Re  =  10,000,  using  Reynolds- averaged  Navier-Stokes 
(RANS)  and  hybrid  RANS/large  eddy  simulation  (LES)  turbulence  modeling  strategies  [Gross 
and  Fasel  2010].  We  also  carried  out  direct  numerical  simulations  (DNS)  of  the  same  square- 
duct  flow  [Gross  and  Fasel  2009a].  Data  from  these  DNS  and  DNS  by  Huser  and  Biringen 
[1993]  served  as  a  reference  for  our  RANS  and  hybrid  RANS/LES  simulations.  For  our  RANS 
and  hybrid  RANS/LES  simulations  we  considered  a  modified  version  of  the  flow  simulation 
methodology  (FSM)  [Speziale  1997,  Fasel  et  al.  2002,  Israel  2005,  Zhang  et  al.  2000]  and  a 
modified  version  of  the  filter-based  RANS  (FBR)  [Johansen  et  al.  2004]  for  which  “backscatter” 
was  introduced  [Gross  and  Fasel  2009a, b].  The  hybrid  models  were  based  on  the  1988  k-w 
model  with  low-Reynolds  number  terms  and  the  1998  k-co  model,  both  with  explicit  algebraic 
stress  model  (EASM)  by  Rumsey  and  Gatski  [2001].  The  FBR  simulations  were  found  to  be  in 
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closer  agreement  with  the  reference  data.  In  particular,  reasonably  accurate  predictions  of  the 
bulk  velocity  and  the  streamwise  velocity  profiles  in  wall  units  were  obtained.  Encouraged  by 
the  channel  flow  results  we  carried  out  hybrid  RANS/LES  simulations  of  the  Stanford  diffuser 
experiment  [Cherry  et  al.  2008]  using  the  FBR  model.  Also,  compared  to  the  experiment,  the 
amount  of  flow  separation  was  underpredicted  in  our  simulations.  However,  the  quality  of  the 
time-averages  was  poor,  since  due  to  the  asymmetric  shape  of  the  diffuser  data  could  only  be 
averaged  in  time.  The  mean  flow  was  found  to  be  strongly  dependent  on  the  underlying 
turbulence  model.  More  simulations  using  longer  time-averages  and  in  particular  also  more 
investigations  using  different  grid  resolutions  are  required  to  further  corroborate  our  hybrid 
results.  In  addition,  our  modified  hybrid  RANS/LES  models  require  further  attention  and 
improvement. 
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